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Abstract 

In  many  settings,  distributed  sensors  provide  dynamic  measurements  over  a  specified 
time  horizon  that  can  be  used  to  reconstruct  information  such  as  parameters,  states 
or  initial  conditions.  This  estimation  task  can  be  posed  formally  as  an  inverse  prob¬ 
lem:  given  a  model  and  a  set  of  measurements,  estimate  the  parameters  of  interest. 
We  consider  the  specihc  problem  of  computing  in  real-time  the  prediction  of  a  con¬ 
tamination  event,  based  on  measurements  obtained  by  mobile  sensors.  The  spread 
of  the  contamination  is  modeled  by  the  convection  diffusion  equation.  A  Bayesian 
approach  to  the  inverse  problem  yields  an  estimate  of  the  probability  density  function 
of  the  initial  contaminant  concentration,  which  can  then  be  propagated  through  the 
forward  model  to  determine  the  predicted  contaminant  held  at  some  future  time  and 
its  associated  uncertainty  distribution.  Sensor  steering  is  effected  by  formulating  and 
solving  an  optimization  problem  that  seeks  the  sensor  locations  that  minimize  the 
uncertainty  in  this  prediction. 

An  important  aspect  of  this  Dynamic  Sensor  Steering  Algorithm  is  the  ability  to 
execute  in  real-time.  We  achieve  this  through  reduced-order  modeling,  which  (for  our 
two-dimensional  examples)  yields  models  that  can  be  solved  two  orders  of  magnitude 
faster  than  the  original  system,  but  only  incur  average  relative  errors  of  magnitude 
0(10“^).  The  methodology  is  demonstrated  on  the  contaminant  transport  problem, 
but  is  applicable  to  a  broad  class  of  problems  where  we  wish  to  observe  certain 
phenomena  whose  location  or  features  are  not  known  a  priori. 

Thesis  Supervisor:  Karen  E.  Willcox 

Title:  Associate  Professor  of  Aeronautics  and  Astronautics 
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Chapter  1 


Introduction 

1.1  Motivation 

In  numerous  fields  of  science  and  engineering  there  is  a  strong  demand  for  precise  ob¬ 
servation  of  large-scale  physical  systems,  modeled  by  partial  differential  equations,  in 
order  to  detect  phenomena  of  interest  that  can  not  be  foreseen  and  to  compute  accu¬ 
rate  predictions  regarding  those  phenomena.  For  example,  let  us  observe  the  amount 
of  contamination  in  the  air  of  the  Boston  city  area  as  measured  by  mobile  sensors 
and  modeled  by  the  convection  diffusion  equation.  Imagine  the  following  scenario: 
due  to  a  terrorist  attack,  a  highly  poisonous  contaminant  is  set  free,  endangering  the 
life  of  many  civilians.  This  situation  calls  for  immediate  results;  counteractions,  like 
evacuations  or  the  deployment  of  the  hre  department,  have  to  take  place  as  soon  as 
possible.  It  becomes  apparent  that  the  time  available  to  come  up  with  an  accurate 
prediction  of  the  ongoing  contamination  is  extremely  limited.  Even  though  the  results 
computed  by  a  high  hdelity,  large-scale  model  are  the  most  accurate  for  this  appli¬ 
cation,  the  long  computational  time  required  by  those  models  (hours  or  even  days) 
makes  them  inappropriate  under  these  circumstances.  Furthermore,  as  the  obtained 
prediction  depends  greatly  on  the  location  of  sensors,  we  need  a  way  to  hnd  the  best 
sensor  locations  for  this  purpose. 

The  discipline  of  model  order  reduction  studies  properties  of  large-scale  dynamical 
systems  and  reduces  their  complexity  while  maintaining  their  input-output  behavior. 
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The  reduced-order  model  is  computationally  much  less  expensive  than  the  full  model 
and  thus  yields  real-time  solutions  for  a  negligible  trade-off  in  accuracy.  The  Bayesian 
approach  to  inverse  problems  is  based  on  observed  data,  which  allow  one  to  estimate 
unknown  model  parameters.  These  two  techniques  will  be  key  ingredients  for  tackling 
the  problem  stated  above.  The  remaining  challenge  is  to  develop  a  methodology  that 
combines  these  two  disciplines.  This  will  enable  accurate  real-time  predictions,  which 
are  optimal  for  the  observed  data  in  systems  that  change  quickly,  and  to  handle 
uncertainty  in  measurements. 


1.2  Literature  Review 

Model  order  reduction  is  a  technique  that  observes  properties  of  large-scale  dynamical 
systems  and  reduces  their  complexity  while  maintaining  their  input-output  behavior. 
The  majority  of  methods  obtains  the  reduced-order  model  by  projecting  the  full 
model  onto  a  basis  that  spans  a  space  of  lower  dimension.  This  reduced  basis  can  be 
computed  using  proper  orthogonal  decomposition  [29,  47],  Krylov-subspace  methods 
[18,  21,  26],  reduced  basis  methods  [42],  balanced  truncation  [36]  or  a  Hessian-based 
model  reduction  approach  [4]. 

The  Bayesian  approach  to  inverse  problems  is  thoroughly  discussed  in  the  liter¬ 
ature,  e.g.  [30,  50]  and  deals  with  estimating  unknown  model  parameters  based  on 
observed  data.  The  solution  of  the  inverse  problem  is  represented  by  a  probability 
density  as  oppose  to  just  a  single  estimate  of  the  model  parameters  as  in  the  determin¬ 
istic  approach.  In  the  Bayesian  approach,  model  as  well  as  measurement  uncertainties 
can  be  incorporated  into  the  solution  very  elegantly. 

In  [34,  35]  an  approach  is  proposed  for  the  description  of  atmospheric  flows  based 
on  proper  orthogonal  decomposition.  In  this  thesis  we  assume  different  velocity  helds 
and  focus  on  the  contaminant  transport  problem.  Combining  the  work  that  was 
presented  in  [34,  35]  and  this  thesis  can  yield  even  more  efficient  and  more  realistic 
results. 

Optimal  sensor  placement  is  an  issue  that  has  been  widely  addressed  in  literature. 
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This  discipline  holds  many  challenges  for  several  different  applications.  The  most 
relevant  papers,  for  the  purpose  of  this  thesis,  shall  be  discussed  now. 

Both,  [6]  and  [44],  address  the  problem  of  placing  sensors  in  water  networks  to 
detect  maliciously  injected  contaminants,  which  has  a  similar  motivation  as  our  prob¬ 
lem.  [6]  formulates  a  linear  mixed-integer  problem  that  minimizes  the  fraction  of  the 
population  at  risk  and  assumes  a  hnite  number  of  sensor  locations.  While  efficiently 
solvable  this  approach,  as  many  others,  only  considers  a  finite  number  of  sensor  lo¬ 
cations.  [44]  focuses  on  hnding  a  layout  of  early  warning  detection  stations  based 
on  an  optimization  problem  comprising  water  quality  conditions  and  extended  pe¬ 
riod  unsteady  hydraulics.  The  main  limitation  of  this  methodology  is  the  real-time 
assumption  and  again  the  hnite  number  of  possible  sensor  locations. 

An  autonomous  model-based  reactive  observing  systems  has  been  developed  in 
[11]  with  a  set  of  static  and  mobile  sensors.  The  focus  lies  on  a  sample  selection 
problem  modeled  by  a  subset  selection  problem  for  regression.  They  assume  the 
physical  phenomena  that  they  wish  to  observe,  do  not  change  temporally  (or  change 
very  slowly).  In  the  contaminant  transport  setting  on  the  other  hand  we  have  to 
be  prepared  for  fast  changing  phenomena.  In  [11]  the  mobile  sensors  are  steered  by 
an  algorithm  based  on  local  linear  regression  as  oppose  to  our  eigenvalue-based  op¬ 
timization  problem.  They  distinguish  between  three  different  kinds  of  sensor  faults 
and  describe  methods  to  detect  those  faults  whereas  we  propose  to  incorporate  mea¬ 
surement  errors  as  uncertainties  when  solving  the  inverse  problem. 

In  [24]  a  sensor  placement  strategy  to  obtain  the  most  effective  visual  sensing 
of  an  area  of  interest  is  proposed  and  solved  by  solving  a  variant  of  the  art-gallery 
problem  [43] .  This  problem  is  NP-hard  and  once  more  the  optimal  solution  is  chosen 
among  a  hnite  number  of  possibilities.  The  issue  of  reducing  uncertainty  over  a 
region  of  interest  has  been  addressed  by  [13]  and  [27].  In  [27]  this  is  approached 
by  an  algorithm  that  maximizes  mutual  information  in  Gaussian  processes  and  [13] 
presents  an  algorithm  for  an  observation  targeting  problem  formulated  as  a  sensor 
selection  problem. 

In  order  to  obtain  the  best  prediction  we  have  to  hnd  the  Bayesian  optimal  exper- 
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imental  design  by  solving  an  optimization  problem  aiming  at  minimizing  uncertainty 
in  the  prediction  based  on  an  approximation  of  the  variance-covariance  matrix.  The 
most  common  design  criteria  focus  on  minimizing  the  trace  of  this  matrix  (this  is 
referred  to  as  A  optimality),  other  approaches  [5,  38]  involve  the  determinant  or 
the  largest  eigenvalue  and  are  known  as  D  and  E  optimal  designs  [2],  The  work 
of  [1,  12,  28]  involves  A  optimal  design.  In  this  thesis  we  focus  on  an  E  optimal 
approach. 

Many  of  the  optimal  sensor  placement  approaches  mentioned  above  compute  an 
optimal  solution  from  a  hnite  number  of  sensor  locations.  In  this  thesis  we  rather 
choose  to  formulate  the  task  as  a  continuous  optimization  problem,  where  the  possible 
sensor  locations  are  assumed  to  be  continuous  and  the  sensors  are  represented  using 
mollihed  delta  functions.  This  formulation  has  the  advantage  that  it  can  be  solved 
efficiently  using  a  gradient-based  algorithm,  enabling  the  possibility  of  real-time  com¬ 
putations. 

1.3  Thesis  Objectives 

The  main  goal  of  this  thesis  is  to  develop  a  Dynamic  Sensor  Steering  Algorithm, 
divided  into  an  offline  and  online  stage,  that  operates  in  real-time  in  order  to  observe 
quickly  changing  phenomena  in  physical  processes.  The  specific  objectives  to  reach 
this  goal  are  the  following: 

1.  As  the  physical  processes  are  modeled  by  partial  differential  equations  and  lead 
to  large-scale  systems,  which  are  too  computationally  expensive  to  be  solved  in 
real-time,  we  propose  an  algorithm  that  builds  a  reduced-order  model  only  once 
in  the  offline  stage  and  repeatedly  runs  through  the  online  stage  using  only  the 
reduced-order  model. 

2.  We  solve  a  Bayesian  formulation  of  the  ill-posed  inverse  problem  to  account  for 
model  and  measurement  uncertainties. 

3.  We  formulate  and  solve  an  optimization  problem  based  on  the  largest  eigen- 
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value  of  the  covariance  matrix,  yielding  optimal  sensor  locations  to  minimize 
uncertainty  in  the  prediction  over  the  output  region  of  interest. 

4.  We  demonstrate  the  newly  proposed  methodology  by  applying  it  to  a  conta¬ 
minant  transport  problem  modeled  by  the  convection  diffusion  equation  and 
extending  the  methodology  to  be  able  to  handle  a  parameterized  velocity  held. 


1.4  Thesis  Outline 

In  Chapter  2  we  develop  the  methodology  for  a  Dynamic  Sensor  Steering  Algorithm 
based  on  a  characterization  of  prediction  uncertainty,  computed  using  a  Bayesian 
formulation  of  the  inverse  problem  and  the  application  of  model  order  rednction. 
Each  component  of  the  approach,  i.e.  the  Bayesian  framework,  the  model  order 
reduction  and  the  optimization  are  discnssed  in  detail. 

In  Chapter  3  the  methodology  developed  in  Chapter  2  is  applied  to  a  contaminant 
transport  problem  modeled  by  the  convection  diffusion  equation.  We  also  include  a 
thorough  comparison  between  fnll  and  rednced-order  model  performance  in  forward 
and  inverse  problem  as  well  as  in  the  optimization. 

We  introdnce  a  parameterized  velocity  held  in  the  Dynamic  Sensor  Steering  Al¬ 
gorithm  setting  in  Chapter  4,  which  enables  a  more  realistic  simulation  of  physical 
processes.  The  posed  challenge  is  to  maintain  the  algorithm’s  computational  real¬ 
time  property,  by  bnilding  reduced-order  models  that  are  not  only  dependent  on  the 
initial  condition  but  also  on  the  velocity  held.  The  resnlts  are  presented  in  Chapter 
5  where  we  employ  a  Navier-Stokes  solver  to  pre-compnte  the  velocity  helds. 

Finally,  Chapter  6  conclndes  the  thesis  with  a  summary  and  snggestions  for  futnre 
work. 
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Chapter  2 


Dynamic  Sensor  Steering 
Algorithm  -  Methodology 


This  chapter  provides  a  description  of  the  methodology  used  in  the  newly  developed 
Dynamic  Sensor  Steering  Algorithm,  which  is  based  on  a  characterization  of  prediction 
uncertainty,  computed  using  a  Bayesian  formulation  of  the  inverse  problem  and  the 
application  of  model  order  reduction.  In  Section  2.1  the  algorithm’s  general  work 
flow,  its  steps  and  its  functionality  are  presented,  followed  by  a  detailed  discussion 
in  Section  2.2  on  how  a  reduced-order  model  is  incorporated  to  obtain  an  accurate 
predictive  capability  that  operates  in  real-time.  We  will  then  analyze  the  Bayesian 
framework  and  the  optimization  problem  employed,  in  Section  2.3  and  Section  2.4 
respectively. 


2.1  Problem  Description 

2.1.1  Introductory  Definitions 

Let  the  two  or  three-dimensional  domain  that  we  wish  to  observe  be  named  D.  In 
order  to  obtain  measurements  of  an  ongoing  physical  process,  Q  mobile  sensors  are 
placed  into  this  domain.  For  the  purpose  of  this  general  discussion  we  do  not  assume 
any  knowledge  about  D.  Constraints  describing  the  sensors’  mobility  have  to  be 


23 


adapted  according  to  each  individual  problem  setting.  Therefore  let  us  collect  the 
corresponding  sensor-related  information  in  the  so-called  sensor  constraints.  One 
approach  to  dehning  these  constraints  is  described  in  Section  2.4.  Let  the  sensor 
locations  be  denoted  as  SI, ...  ,Sq  for  i=l,2,3. . . ,  where  the  superscript  refers  to  a 
time  index  and  the  subscript  corresponds  to  the  number  of  the  sensor.  Furthermore 
let  Tsteer  deuote  the  available  steering  time  for  the  mobile  sensors.  We  do  not  assume 
anything  about  the  initial  sensor  placement  except  that  they  have  to  be  within  our 
domain  O.  Within  this  domain  there  are  some  specihc  regions  that  are  of  greater 
interest  than  others.  It  may  be  temporary  interest  because  sensors  detect  some  curious 
unexpected  phenomenon  or  permanent  interest  due  to  additional  knowledge  on  the 
structure  of  hi,  e.g.  residential  areas.  The  union  of  those  regions  will  be  referred  to  as 
the  output  region  of  interest  and  be  denoted  as  flj.  Let  g{-)  be  the  forward  operator 
used  to  model  the  natural  processes  we  want  to  observe  e.g.  the  convection-diffusion 
equations,  let  m  be  the  input  vector  for  the  model  g,  and  let  d  denote  the  vector 
containing  the  outputs  of  interest.  Then  this  can  be  written  as  follows 

d  =  g{m).  (2.1) 

In  the  case  that  g  is  a  linear  operator,  equation  (2.1)  can  be  written  as: 

d  =  Gm.  (2.2) 

In  Section  2.2,  we  will  present  a  model  reduction  framework.  For  now,  we  define  the 
following  notation.  Let  fl'r(-)  denote  the  reduced-order  model  for  g{-)  and  Gr  be  the 
reduced-order  model  for  G,  then  we  obtain  the  two  reduced  systems  (2.3)-(2.4)  that 
correspond  to  (2.1)  and  (2.2),  respectively.  The  goal  of  model  order  reduction  is  to 
determine  a  reduced  model  so  that  d  k,  d^  while  achieving  a  signihcant  reduction  in 
the  computational  cost. 


dr 

=  grim), 

(2.3) 

dr 

=  GrTtn 

(2.4) 
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2.1.2  Contaminant  Transport  Example 

To  demonstrate  the  results  of  the  Dynamic  Sensor  Steering  Algorithm  we  consider  a 
two-dimensional  contaminant  transport  problem,  modeled  by  the  convection-diffusion 
equation, 


-|-u-Vu  — = 

0 

in  D  X  (0, f/). 

(2.5) 

u  = 

0 

on  Tt?  X  (0,t/), 

(2.6) 

(9u 

dn 

0 

on  Ttv  X  (0,t/), 

(2.7) 

u  = 

Uo 

in  D  for  t  =  0, 

(2.8) 

where  u  denotes  the  contaminant  concentration,  v  denotes  the  velocity  held,  tj  de¬ 
notes  the  observed  time  horizon  and  k  is  the  diffusivity  constant.  We  impose  homo¬ 
geneous  Dirichlet  boundary  conditions  on  the  Dirichlet  boundary  T x  (0,t/)  and 
Neumann  boundary  conditions  on  the  Neumann  boundary  T tv  x  (0,t/).  We  perform 
the  time  discretization  using  the  Backward  Euler  method  and  a  Streamline  Upwind 
Petrov-Galerkin  [8]  hnite-element  method  using  triangular  elements  for  the  discretiza¬ 
tion  in  space  to  obtain  the  discrete-time  system 

Di  u{k  -|-  1)  =  D2  u{k),  /c  =  0, 1, . . . ,  T  —  1  (2.9) 

d{k)  =  C  u{k),  k  =  0,l,...,T  (2-10) 

w(0)  =  uo,  (2.11) 

where  u{k)  G  represents  the  state  at  time  tk,  d{k)  G  is  the  output  at  time  tk 
and  the  initial  condition  uq  is  the  state  at  time  t  =  0  and  the  matrices  Di  G 
D2  G  C  G 

2.1.3  The  Algorithm 

The  Dynamic  Sensor  Steering  Algorithm  focuses  on  computing  the  best-possible  pre¬ 
diction  of  an  ongoing  physical  process  in  the  output  region  of  interest  D/  in  real-time 
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based  on  sensor  measurements.  In  order  to  improve  the  amount  of  information  ob¬ 
tained  from  the  measurements  and  thereby  increase  the  accuracy  of  the  prediction,  we 
solve  an  optimization  problem  to  determine  new  sensor  locations  that  lead  to  mini¬ 
mized  uncertainty  in  the  prediction.  The  algorithm  is  made  up  of  two  separate  stages, 
the  offline  stage,  which  we  execute  only  once,  and  the  online  stage,  which  consists  of  a 
measure-predict-optimize-steer  cycle.  This  online  cycle  is  repeatedly  executed  to  im¬ 
prove  the  accuracy  of  the  obtained  prediction,  while  incorporating  new  measurement 
information  as  time  progresses.  The  formation  of  a  reduced-order  model  in  the  offline 
stage  enables  the  successful  application  of  the  Dynamic  Sensor  Steering  Algorithm  to 
problems  that  require  real-time  processing.  The  online  stage  consists  of  hve  major 
steps  computing  a  current  prediction,  solving  an  optimization  problem  on  the  sensor 
locations  and  steering  each  mobile  sensor  to  its  new  and  improved  destination.  Note 
that  for  this  section  the  input  vector  m  consists  of  an  initial  condition  denoted  by  uo 
only. 


Offline  Stage 


step  Oa.  Build  reduced-order  model  gr(uo)  of  physical  system  g(uo). 


Step  Ob.  Initialize  sensor  locations  Si°,.--,Sq°. 


Figure  2-1:  Offline  stage  of  the  Dynamic  Sensor  Steering  Algorithm. 

Offline  Stage  of  the  Dynamic  Sensor  Steering  Algorithm 

Figure  2-1  shows  the  offline  stage  of  the  algorithm.  The  details  of  each  step  are 
discussed  in  the  following: 

Step  Oa.  Build  reduced-order  model  gr{-)  of  the  physical  system  g{-)  and  continue 
working  with  reduced  system  throughout  the  Dynamic  Sensor  Steering  Algorithm. 
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Note  that  some  well-known  techniques  to  construct  reduced-order  models  and  the 
method  applied  in  this  thesis  for  the  application  to  a  contaminant  transport  problem 
are  discussed  in  Section  2.2. 

Step  Ob.  Place  the  Q  mobile  sensors  at  locations  S^, . . . ,  Sq  within  hi  and  acti¬ 
vate  them.  We  note  that  specihcation  of  the  initial  sensor  locations  is  an  important 
question;  however,  in  the  setting  of  this  thesis  we  focus  on  the  dynamic  element,  i.e. 
mobile  sensors  that  are  moved  to  optimal  locations  so  as  to  minimize  uncertainty  in 
the  prediction. 

Online  Stage  of  the  Dynamic  Sensor  Steering  Algorithm 

Figure  2-2  shows  the  online  stage  of  the  algorithm.  The  details  of  each  step  are 
discussed  in  the  following: 

Step  1.  The  Q  sensors  create  measurements  of  the  environment  at  their  current 
locations  . . . ,  in  cyele  i.  Let  (P  =  . . . ,  (Pq)  be  the  obtained  data  at  the 

Q  different  sensor  locations  in  cycle  i. 

Step  2.  Based  on  the  data  d*  =  {d\, . . .  jdq)  and  the  known  forward  model  fl'r(') 
{Gr  in  the  linear  case),  we  solve  the  inverse  problem  using  a  Bayesian  approach  and 
compute  the  posterior  probability  density  Please  refer  to  Section  2.3  for 

further  details  on  the  Bayesian  framework. 

Step  3.  In  many  applications  some  regions  12/  within  the  observed  domain  hi  need 
to  be  watched  more  thoroughly  than  others  at  any  time,  meaning  that  the  physical 
process  that  we  want  to  predict  interests  us  especially  in  those  regions.  12/  can  be 
hxed,  e.g.  representing  an  urban  region,  or  could  change  from  cycle  to  cycle,  e.g.  if 
we  wish  to  track  regions  of  high  concentration. 

Step  4-  We  optimize  the  sensor  locations  within  12  such  that  the  uncertainty 
in  the  prediction  for  the  output  regions  of  interest  12/  is  minimized.  For  details  on 
how  the  optimization  problem  is  set  up  and  solved  please  refer  to  Section  2.4.  As 
input  parameters  we  use  the  current  set  of  sensor  locations  a  set  of 

constraints  specifying  allowable  sensor  movements,  knowledge  about  the  domain  12 
and  12/,  and  the  model  -  either  g  oi  Qr-  As  the  output  of  this  optimization  problem 


27 


we  obtain  the  optimal  sensor  locations  S{, ...  ,Sq,  whose  quality  will  be  reflected  in 
Step  2  of  the  next  cycle  of  the  online  stage,  where  a  new  prediction  will  be  produced. 

Step  5.  The  last  step  in  the  algorithm  steers  the  mobile  sensors  from  their  cur¬ 
rent  locations  ,  S'q^  to  the  optimized  locations  5^, . . . ,  S'g  within  Tgteer  and 

return  to  Step  1.  Note  that  the  choice  of  Tgteer  is  connected  to  the  above  mentioned 
constraints  on  sensor  movements.  Depending  on  the  application,  there  may  be  a 
trade-off  between  extending  Tgteer,  thereby  increasing  the  amount  of  reachable  sensor 
locations  and  limiting  the  number  of  cycles  that  can  be  run  in  the  time  available, 
or  rather  focusing  on  gaining  as  many  predictions  as  possible  and  restricting  sensor 
movements. 


2.2  Model  Reduction  Framework 

2.2.1  Significance  for  Dynamic  Sensor  Steering  Algorithm 

In  general,  model  order  reduction  is  applied  to  highly  complex  large-scale  dynamical 
systems  to  obtain  a  reduced-order  model  (ROM)  whose  dimensions  are  considerably 
smaller,  thus  leading  to  faster  computational  time  while  the  input-output  behavior 
remains  approximately  the  same.  In  the  majority  of  cases  model  reduction  techniques 
project  the  large-scale  systems  onto  a  basis  in  a  space  of  reduced  dimension.  Some 
well-known  methods  to  obtain  this  reduced  basis  are  proper  orthogonal  decomposition 
(POD)  [14,  29,  47],  approximate  balanced  truncation,  Krylov-subspace  methods  [18, 
21,  26]  and  Hessian-based  model  reduction  [4]  for  initial  value  problems. 

Let  us  emphasize  once  more  the  significance  of  model  order  reduction  for  realistic 
applications  of  the  Dynamic  Sensor  Steering  Algorithm.  In  order  to  run  the  online 
stage  in  real  time,  a  ROM  of  the  physical  model  has  to  be  computed  offline  making 
sure  computations  remain  accurate  but  are  considerably  faster.  Figure  2-3  depicts 
the  trade-off  between  full  and  reduced-order  model. 

After  construction  in  the  offline  stage  the  reduced-order  model  is  used  at  two 
points  within  the  online  stage:  in  Step  2  when  solving  the  inverse  problem  and  in 


Online  Stage 
Cycle  i 


Figure  2-2:  One  cycle  of  the  online  stage  of  the  Dynamic  Sensor  Steering  Algorithm. 
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Step  4  when  solving  the  optimization  problem.  For  a  detailed  comparison  between  full 
and  reduced-order  model  performance  in  a  two-dimensional  contaminant  transport 
problem  please  refer  to  Chapter  3. 


2.2.2  Model  Reduction  Methodology 

In  this  section  let  us  focus  on  how  to  build  the  reduced-order  models  intended  for  this 
thesis.  Therefore  let  us  restate  the  linear  discrete-time  system  describing  an  initial 
value  problem  from  Section  2.1.2  as 


Di  u{k  +  1) 

=  D2  u{k), 

A;  =  0,1,. 

..,T-  1 

(2.12) 

d{k) 

=  C  u{k), 

A;  =  0,1,.. 

■,T 

(2.13) 

«(0) 

=  UO: 

(2.14) 

where  u{k)  G  represents  the  state  at  a  time  tk  ior  k  =  0, . . .  ,T,  d{k)  G  is  the 
output  at  time  t  and  the  initial  condition  uq  is  the  state  at  time  t  =  0.  We  obtain 
Di,D2  G  and  C  G  from  spatial  and  temporal  discretization  methods. 

Then  the  system  from  (2.12)  can  be  rewritten  in  matrix  form 

Au  =  F  uo,  (2.15) 

d  =  C  u,  (2.16) 


where 


(2,17) 
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S'~Sr' 
m'  ~  nir’ 


Real-time 

computations 


_ * _ 

Output: 

m/,  S>(Si;,...,Sq;) 


Figure  2-3:  Trade-off  between  full  and  reduced-order  model  in  one  online  cycle. 


31 


and  the  matrices’  dimensions  are  given  by  A  =  F  =  and 

C  =  and  where  u  and  d  are  dehned  as 


u(0) 

d{0) 

«.(1) 

,  d  = 

d{l) 

u(T)_ 

d{T)_ 

(2.18) 


As  N  is  usually  large  we  need  to  create  a  reduced-order  model  by  projecting  the 
large-scale  system  onto  a  basis  in  a  space  of  considerably  reduced  dimension  n. 
Let  V  G  denote  the  matrix  containing  the  n  orthonormal  basis  vectors  and 

Ur(k)  G  M”  denote  the  reduced-order  state,  then  Vur{k)  should  provide  us  with  a 
good  approximation  for  the  state  u{k)  for  /c  =  0, 1, . . . ,  T.  Then  we  apply  a  Galerkin 
projection  method  to  the  full  system  by  projecting  onto  the  space  spanned  by  V  and 
hnally  obtain  the  reduced-order  model  of  the  discrete-time  system  as 


Uy{k  +  1) 

=  D2,y  Uy{k),  k  =  D,l, . 

..,T-1 

(2.19) 

dy{k) 

=  Cy  Uy{k),  /c  =  0,  1,  .  . 

■,T 

(2.20) 

Mr(0) 

=  uo. 

(2.21) 

where  =  V^DiV  for  i  =  1,2  and  Cr  =  CV .  Note  that  this  can  also  be  rewritten 
in  matrix  form 


Ay.  Uy  =  Fy  UQ  ,  (2.22) 

dy  =  Cy  Uy,  (2.23) 

where  Ay,  Cy,  Uy  and  dy  are  defined  exactly  like  A,  C,  u  and  d  in  (2.17)  &  (2.18) 
except  for  the  fact  that  each  inner  component  of  the  vectors  or  matrices  has  a  subscript 
r  now  and  thus  we  only  have  to  exchange  iV  to  n  to  get  the  dimensions  of  the  reduced 
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system.  Only  G  is  defined  slightly  differently,  as 


1/^ 

0 


0 


(2.24) 


We  use  the  Hessian-based  model  reduction  approach  in  [4]  where  s  seed  initial  con¬ 
ditions  are  computed  as  the  dominant  eigenvectors  of  the  system’s  Hessian  dehned 
as 

H  =  {CA-^Ff{CA-^F).  (2.25) 

Then  state  trajectories  at  time  steps  k  =  0, 1,2, ...  ,T  are  generated  for  each  seed 
initial  condition  and  stored  in  the  snapshot  matrix  X  G  where  M  =  s{T  +  1). 

The  technique  used  to  efficiently  compute  the  basis  V  from  all  the  state  trajectories 
is  proper  orthogonal  decomposition  (POD).  There  are  several  ways  to  compute  the 
POD  [31]  but  for  our  means  we  will  concentrate  on  a  method  using  the  singular  value 
decomposition  of  the  snapshot  matrix  X.  Let  the  singular  value  decomposition  for 
the  matrix  of  snapshots  be  dehned  as 


X  =  USV'^, 

then  the  following  holds 

=  US^U^, 

x^x  =  vs^v^. 


(2.26) 


(2.27) 

(2.28) 


Due  to  the  above  equations,  the  singular  values  of  X  are  the  square  roots  of  the 
eigenvalues  of  X^X  or  XX^  and  moreover  the  left  and  right  singular  vectors  of  X 
are  in  fact  the  eigenvectors  of  XX^  and  X^X.  The  proper  orthogonal  modes  are 
the  left  singular  vectors  of  X  and  the  proper  orthogonal  values  are  dehned  as  the 
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Observed  data 

d 


Forward  Problem 


Solve  for  d 


by  using  g(.),  m 


Inverse  Problem 


. I . 

Solve  for  m 

III 

by  using  g(.),  d 


Figure  2-4:  Comparison  between  inverse  and  forward  problem. 

corresponding  singular  values.  As  we  started  out  to  gain  a  basis  containing  n  vectors, 
we  just  have  to  choose  the  dominant  n  left  singular  values  of  X  to  obtain  the  basis. 

2.3  Bayesian  Inverse  Approach 

Throughout  this  section  we  are  following  the  notation  introduced  by  [22,  50].  and 
thus  let  g{-)  denote  the  physical  system  we  wish  to  observe,  let  m  denote  the  model 
parameters  and  let  d  be  the  observable  parameters,  such  that  the  theoretical  forward 
problem  can  be  written  as 

d  =  g{m).  (2.29) 

In  each  online  cycle  of  the  Dynamic  Sensor  Steering  Algorithm  we  obtain  measure¬ 
ments  d*  that  are  related  to  model  parameters  m*  via  the  physical  model  g{-)  as 
described  in  (2.29).  In  order  to  predict  m*  using  the  observed  data  d*  and  the  for¬ 
ward  operator  g{-)  an  inverse  problem  needs  to  be  solved.  While  the  forward  problem 
maps  from  the  model  parameters  or  the  “cause”  to  the  data,  the  inverse  problem  solves 
for  the  model  parameters  using  the  data  and  the  knowledge  of  the  physical  model, 
see  Figure  2-4.  In  practice  the  theoretical  results  for  d  obtained  by  using  (2.29)  will 
never  exactly  be  the  same  as  the  actual  observations  that  are  made  during  the  on- 
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line  stage  of  the  Dynamic  Sensor  Steering  Algorithm  due  to  model  uncertainties  and 
measurement  errors.  Thus,  when  reconstructing  m  based  on  g{-)  and  d,  it  is  hard  to 
quantify  whether  the  obtained  solution  m  is  correct.  What  we  are  rather  interested 
in  is  a  way  of  saying  that  some  solutions  are  in  fact  more  likely  than  others.  In  order 
to  accurately  incorporate  the  model  and  measurement  uncertainties  and  the  likeliness 
of  a  solution  m  we  introduce  probability  density  functions  to  the  inverse  problem, 
leading  us  the  Bayesian  inverse  approach  where  we  compare  theoretical  predictions 
to  observations.  One  basic  tool  needed  in  this  framework  is  Bayes’  Theorem  relating 
cause  to  effect. 


2.3.1  Bayes  Theorem 

Let  A  and  B  be  two  events,  then  the  joint  probability  can  be  written  as 

P(A,  B)  =  P{A  I  B)  P{B)  =  P{B  I  A)  P{A),  (2.30) 

where  P{A  \  B)  is  the  conditional  probability  of  A  in  case  event  B  has  already 
happened.  From  this  we  immediately  obtain  Bayes’  Theorem 

I  p, .  „3i) 

Let  us  now  incorporate  Bayes’  Theorem  into  the  Dynamic  Sensor  Steering  setting  in 
the  following  section. 

2.3.2  Solution  of  the  Inverse  Problem 

Let  M.  denote  the  model  space  manifold  and  T>  denote  the  data  space  manifold. 
Furthermore  let  and  fJ,D{d)  denote  the  homogeneous  probability  density  of 

the  model  space  manifold  and  the  data  space  manifold  respectively,  then  their  joint 
homogeneous  probability  density  can  be  written  as  follows 

fx{d,m)  =  laoid)  iaM{m).  (2.32) 
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Let  6{d\m)  be  the  conditional  probability  density  that  provides  knowledge  about  d 
when  we  have  m.  This  is  determined  by  the  forward  operator  and  hence  called  the 
likelihood  function.  The  joint  probability  density  correlating  m  and  d  is  dehned  as 

Q{d,m)  =  9{d\m)  jj,M{'m).  (2.33) 


Due  to  the  fact  that  there  are  measurement  and  model  uncertainties  let  us  introduce 
the  joint  prior  probability  density  p{d,m).  When  solving  the  inverse  problem  we 
are  interested  in  the  information  which  will  be  provided  by  the  posterior  density.  In 
particular,  we  are  interested  on  the  posterior  information  on  the  model  parameters 
m.  From  Bayes’  Theorem  we  obtain  the  joint  posterior  probability  density  ^{d,  m) 
as  follows: 


a{d,  m)  =  q 


Q{d,  m)  p{d,  m) 
p,{d,  m) 


(2.34) 


where  g  is  a  normalizing  constant.  From  (2.34)  we  obtain  the  two  marginal  probability 
densities  and  aoid)  as  follows: 


aM{rn)  = 

/  (T{d,m] 
Jm 

1  dm 

(2.35) 

CToid)  = 

/  a{d,  m) 

Jv 

dd 

(2.36) 

(Tm(^)  contains  the  posterior  information  on  the  model  parameters  given  the  ob¬ 
served  data  and  is  therefore  the  solution  to  the  inverse  problem  in  the  Bayesian 
setting.  Due  to  the  fact  that  the  number  of  model  parameters  may  be  large  it  might 
be  difficult  to  be  able  to  extract  information  from  (TA^(m).  Thus  in  order  to  analyze 
this  posterior  probability  density  in  more  detail  we  want  to  observe  marginal  densi¬ 
ties.  Let  m  =  (mi, . . .  ,mp)  where  P  denotes  the  number  of  model  parameters  and 
let  M.  =  M-i  X  ...  X  M.p  he  the  Cartesian  product  of  P  model  manifolds,  then  the 
posterior  marginal  probability  density  for  m*  for  an  1  <  i  <  P  can  be  obtained  by 


a. 


rrii 


.nii  = 


drui . . . 


drui-i 


dm 


i+l 


dmp  (Tiw(m).  (2.37) 


'Ml 


'M^-i 


'Mi+i 


'Ml 
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2.3.3  Gaussian  Case 


This  section  is  dedicated  to  restating  the  posterior  probability  density  defined  in 
Section  2.3.2  in  case  uncertainties  are  Gaussian.  Let  us  assume  theoretical  model¬ 
ing  uncertainties,  described  by  the  covariance  matrix  are  Gaussian.  Then  the 
conditional  probability  6{d\m)  relating  m  and  d  becomes 

6{d\m)  ~  exp{—^{d  —  g{m))'^C^^{d  —  g{m))).  (2.38) 

Now  let  Cd  denote  the  measurement  uncertainty  matrix  and  Cm  the  prior  uncer¬ 
tainty  matrix,  then  let  Cd  =  Cd  +  Ct  denote  the  matrix  combining  modeling  and 
measurement  uncertainties.  The  resulting  posterior  probability  density  is  as  follows: 

(JM{m)  ~  exp{-]^{g{m)  -  -  d)  -  ^(m  -  rUpriorY {m  -  niprior)) 

(2.39) 

If  the  forward  problem  is  linear  the  term  g{m)  will  be  substituted  by  Gm  and  the 
posterior  model  covariance  Cm  is  the  inverse  of  the  Hessian  Cm  =  {C^C^^G  + 

Cm")-^- 

2.4  Optimization  Problem 

As  already  mentioned  in  Section  2.2.2,  a  discrete-time  system  can  be  rewritten  in 
matrix  form  and  we  obtain  a  system  such  as  that  stated  in  (2.15)-(2.16).  Let  us 
consider  the  linear  case  where  we  can  state  the  forward  problem  with  forward  operator 
Gi  =  CiA~^F  as  follows 


A  u  =  F  uq,  (2.40) 

d  =  CiU,  (2.41) 

where  (as  before)  u  is  the  space-time  state,  Uq  is  the  initial  condition,  and  d  are  the 
observables  in  space-time.  The  matrix  Ci{S)  defines  the  current  sensor  locations  S  in 
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the  domain  hi.  Due  to  simplicity  of  notation  we  will  be  writing  Ci  instead  of  Gj(S'). 
Similarly  the  prediction  problem  with  prediction  operator  Gp  =  CpA~^F  is  given  by 

A  u  =  F  uo,  (2.42) 

p  =  CpU,  (2.43) 

where  p  are  the  prediction  outputs  of  interest,  dehned  by  the  matrix  Cp.  The  pre¬ 
diction  and  the  inversion  problem  are  formulated  separately  because  in  general  they 
are  not  the  same,  e.g.  sensor  locations  and  locations  of  interest  might  not  coincide  or 
the  problems  have  a  differing  time  horizon. 

Let  denote  the  uncertainty  matrix  combining  model  and  measurement  un¬ 
certainties  and  denotes  the  priori  uncertainty  matrix.  Then  the  expression  for 
the  posterior  of  the  initial  condition  covariance  is 

Cm  =  Cm{S)  =  (F^A-'^Cf  +  C^i)-^(2.44) 

Assuming  is  the  identity  and  0,  which  signihes  that  there  is  no  prior 

information,  we  compute  the  simplihed  posterior  model  covariance  as 

Cm  =  Cm{S)  =  {F^A-'^Cj  C,A~^F  +  (31)-^  =  {GJG,  +  (3I)-\  (2.45) 

Note  that  j3I,  with  (3  small,  is  merely  a  regularization  term  needed  for  inversion. 
The  goal  addressed  in  the  optimization  problem  is  to  minimize  the  uncertainty  in 
the  prediction,  thus  the  posterior  covariance  of  the  prediction  held  will  be  observed, 
computed  as 


Cp  =  Cp{S)  =  GpCUS)Gl.  (2.46) 

With  p,  the  prediction  outputs  of  interest,  being  known,  the  posterior  covariance  of 
the  prediction  held  depends  on  the  set  of  sensor  locations  S.  Thus  we  have  to  move 
the  sensors  to  locations  such  that  the  uncertainty  is  minimized.  This  can  be  achieved 
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by  minimizing  the  spectral  norm  of  Cp,  which  in  turn  is  equivalent  to  minimizing 
the  maximum  eigenvalue  of  Cp.  Therefore  the  dynamic  sensor  steering  optimization 
problem  reads 

min  Amax((^p(*§))-  (2.47) 

s&n 

As  we  are  working  with  mobile  sensors,  we  need  to  impose  constraints  representing 
allowable  sensor  motions.  Those  sensor  constraints  are  incorporated  into  the  opti¬ 
mization  problem,  meaning  that  prediction  uncertainty  will  be  minimized  taking  into 
consideration  only  the  sensor  locations  that  can  be  reached  within  Tgteer  and  that 
satisfy  the  imposed  sensor  constraints.  Thus  the  constrained  optimization  problem 
reads 

min  \ma^{Cp{S))  (2.48) 

s&n 

s.t.  S  satisfies  sensor  constraints 

A  realistic  and  for  this  setting  very  useful  approach  incorporates  the  previously  dis¬ 
cussed  steering  time  Tgteer  in  which  the  sensors  move  to  the  recently  computed  optimal 
locations  that  minimizes  uncertainty  in  the  prediction.  Obviously  we  do  not  want 
Tsteer  to  be  too  long  as  the  Dynamic  Sensor  Steering  Algorithm  aims  at  achieving 
real-time  measuring-predicting-optimizing-steering,  where  time  lost  in  steering  the 
sensors  is  needed  badly  during  the  optimization  as  well  as  during  the  computation 
of  the  prediction.  Thus  we  assume  that  a  sensor  can  only  reach  locations  that  are 
within  a  certain  radius  R  within  the  available  steering  time,  which  can  be  formulated 
as 

min  Amax(Cp(A))  (2.49) 

s&n 

s.t.  - Sjf  +  (S]  - Vj  =  i,...,g, 

where  S  =  {Si, . . . ,  Sq)  is  the  set  of  current  sensor  locations  and  each  sensor  location 
consists  of  an  x  and  y  coordinate,  i.e.  Sj  =  (Sj,  Aj). 

We  solve  the  optimization  problem  using  a  gradient-based  method,  i.e.  a  se¬ 
quential  quadratic  programming  (SQP)  method.  At  each  iteration  in  this  method 
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a  quadratic  programming  subproblem  is  solved,  the  estimate  of  the  Hessian  of  the 
Lagrangian  is  updated  by  the  Broyden-Fletcher-Goldfarb-Shanno  (BFGS)  formula 
[9,  19,  23,  46], 

2.5  Implementation  Remarks 

In  large-scale  problems  there  arise  various  computational  challenges  for  the  newly 
introduced  Dynamic  Sensor  Steering  Algorithm.  First,  the  resulting  model  matrices 
may  have  large  dimensions.  Moreover  as  the  covariance  matrix  is  the  inverse  of  the 
Hessian,  we  need  to  compute  this  inverse  and  its  eigenvalue  spectra  as  well.  Therefore 
it  is  important  to  avoid  explicit  formation  of  those  matrices  and  apply  matrix-free 
methods  only  requiring  the  generation  of  matrix-vector  products. 

For  large-scale  problems  with  high-dimensional  parameter  spaces,  computing  co- 
variance  matrices  can  become  a  very  demanding  problem  because  explicit  formation 
could  already  exceed  available  storage.  Furthermore  to  compute  the  covariance  ma¬ 
trix  the  Hessian  needs  to  be  inverted  (see  Section  2.4)  which  could  be  too  costly. 
Therefore  we  aim  at  approximating  Cm  by  using  both  a  truncated  spectral  decom¬ 
position  i.e.  a  Lanczos  method  and  the  Sherman-Morrison- Woodbury  formula  for 
inversion. 

Let  V KV'^  be  the  eigenvalue  decomposition  of  GfGj, where  A  is  the  diagonal  matrix 
containing  the  eigenvalues  and  V  contains  the  corresponding  eigenvectors.  It  often 
occurs  that  the  eigenvalue  spectrum  of  Gf  Gi  decays  rapidly,  enabling  the  approxima¬ 
tion  through  a  truncated  spectrum,  meaning  that  we  can  approximate  V NV'^  based  on 
the  few  dominant  eigenvalues  and  vectors  contained  in  W,  such  that  V NV'^  ~  VrAV^. 
Applying  the  Sherman- Morrison- Woodbury  formula  [52]  now  we  obtain  the  following 
approximation 

Cm  =  iGjG,+pi)~^  =  {VAV^+pi)-^  ^  (KAH7+/3J)-i  =  (2.50) 

where  D  is  a  diagonal  matrix  with  Da  = 
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A  second  challenge  is  related  to  representation  of  sensor  locations.  In  order  to 
make  the  problem  differentiable,  the  sensor  locations  are  represented  by  mollihed 
delta  functions,  meaning  the  sensor  locations  are  modeled  by  Gaussian  functions 
centered  at  the  according  sensor  location  where  the  variance  of  the  Gaussian  needs 
to  be  chosen  appropriately  in  conjunction  with  the  local  grid  size,  thereby  enabling 
the  sensor  locations  to  be  anywhere  in  the  domain  G  (not  just  at  grid  points). 
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Chapter  3 


Dynamic  Sensor  Steering 
Algorithm  -  Application  and 
Results 


In  Chapter  3  we  apply  the  methodology  presented  in  Chapter  2  to  a  two-dimensional 
contaminant  transport  problem  that  is  modeled  by  the  convection-diffusion  equa¬ 
tion.  This  problem  is  introduced  in  Section  2.1.2.  The  following  sections  present 
results  that  include  a  comparison  between  the  full  and  the  reduced-order  model,  and 
a  demonstration  of  the  use  of  model  reduction  for  uncertainty  quantihcation. 


3.1  Description  of  Rectangular  Domain 


Let  us  consider  a  specihc  example,  where  the  dimensions  of  the  domain  12  are  described 
by 


0  <  a;  <  1 

12  =<  (3.1) 

0  <  y  <  0.4 


and  the  computational  domain  has  spatial  mesh  size  of  iV  =  4005.  Figure  3-1  shows 
the  computational  mesh  with  triangular  elements  as  well  as  the  dimensions  of  domain 
12.  At  the  inflow  boundary  characterized  by  T^)  =  {x  =  0, 0  <  y  <  0.4}  we  impose 
homogeneous  Dirichlet  boundary  conditions.  On  all  the  remaining  boundaries  of  12 
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we  impose  homogeneous  Neumann  boundary  conditions.  Throughout  this  chapter  we 
assume  a  diffusivity  n  =  0.01,  At  =  0.005,  forty  time  steps  and  a  constant  velocity 
held  vi,  dehned  as  follows 


vi{x,y) 


1 

0 


for  (x,  y)  G  hi. 


The  output  region  of  interest,  hi/,  is  dehned  as  follows: 


(3.2) 


Vtj 
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Figure  3-1:  Triangular  mesh  of  domain  hi. 


3.2  Reduced- Order  Model  Performance  in  Forward 
Problem 


We  choose  to  demonstrate  the  results  using  a  superposition  of  three  Gaussian  func¬ 
tions  as  the  initial  contaminant  concentration  depicted  in  Figure  3-2.  Each  Gaussian 
function  is  dehned  as 


Uo{.x,y) 


(3.4) 
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where  a  corresponds  to  the  standard  deviation  and  (xc,  Vc)  represents  the  center  of 
the  Gaussian. 


0.4 


Contaminant  Concentration  at  t=0 


0.3 

.!2 

™  0.2 

>. 


0.1 


0^ - 

0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1 

x-axis 


0.01  0.02  0.03  0.04  0.05  0.06  0.07  0.08  0.09 


Figure  3-2:  This  hgure  depicts  the  sample  test  initial  condition  used  to  compare 
reduced  model  to  full  model  in  the  forward  problem. 

In  Figure  3-3  we  compare  reduced  model  outputs  to  full-scale  outputs  at  two 
sensor  locations  Si  =  (0.15,0.25)  and  S2  =  (0.45,0.2).  We  dehne  the  error  e:  due  to 
a  particular  initial  condition  uo  as 

e:  =  ||d-d^||2  =  \\  {C A~^F  -  CrA~^Fr)  uq  \\2,  (3.5) 

where  d  and  dr  correspond  to  the  full  model  output  and  the  reduced-order  model 
output  respectively.  Table  3.1  displays  various  reduced-order  model  properties,  where 
£  denotes  the  error  compared  to  the  full  model,  c  —  time  denotes  the  computational 
time  to  obtain  the  outputs  in  seconds,  n  denotes  the  size  of  the  reduced-order  model, 
and  the  signihcance  of  A  and  /2  are  as  follows.  When  creating  a  snapshot  matrix  using 
the  Hessian-based  method  described  in  [4]  we  compute  s  eigenvectors  corresponding 
to  the  s  largest  eigenvalues  of  the  Hessian.  Thus  we  specify  A  so  that  eigenvector 
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i,  which  has  eigenvalue  Aj,  is  in  the  pool  of  initial  conditions  if  ^  >  A,  where  Ai  is 
the  largest  eigenvalue  of  the  Hessian.  Then  we  apply  POD  to  the  snapshot  matrix 
to  obtain  a  basis  of  size  n  for  the  reduced-order  model.  The  number  n  is  specihed 
by  relating  the  POD  eigenvalues  pi, . . . ,  to  the  largest  one  pi.  Thus  the  POD 

basis  vector  corresponding  to  eigenvalue  pj  is  included  in  the  basis  V  ii  ^  >  Jl.  Note 
that  the  size  of  the  reduced-order  model  used  in  our  results  will  vary  according  to  the 
choice  of  A  and  Jl. 

Running  the  forward  problem  using  the  full  model  takes  5.5  seconds  on  a  desktop 
computer  with  Intel  Core  2  Duo  processor  and  2GB  RAM.  Note  that  all  computational 
results  in  this  thesis  were  computed  using  the  same  desktop  computer.  From  Table  3.1 
we  see  that  the  largest  reduced-order  model  of  size  n  =  212  is  about  ten  times  faster 
than  the  full  model  and  the  smallest  reduced-order  model  of  size  n  =  38  computes 
results  360  times  faster.  The  performance  with  respect  to  error  and  computational 
time  of  each  reduced-order  model  can  be  seen  in  Figure  3-4. 
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118 

0.1 

10-6 

5 

1.9211e-^ 

0.1480 

121 

0.01 

o 

1 

6 

1.9087e-^ 

0.1950 

135 

0.001 

10-^ 

7 

1.1382e-^ 

0.4340 

188 

0.01 

10-6 

8 

1.0232e-^ 

0.5690 

212 

0.001 

10-6 

Table  3.1:  Properties  of  various  reduced-order  models  of  a  full-scale  system  with  size 
N  =  4005  and  two  output  sensors  in  the  forward  problem.  The  error  e  is  dehned  in 
(3.5)  and  the  initial  condition  used  is  depicted  in  Figure  3-2. 


3.3  Solution  of  the  Inverse  Problem 

In  the  contaminant  transport  problem  introduced  in  (2.5)-(2.8)  the  vector  of  model 
parameters  m  corresponds  to  the  initial  condition  uq,  and  the  output  or  the  data 
parameters  are  d  =  {di, ,  dq),  where  dj  ior  1  <  j  <  Q  corresponds  to  the  conta- 


46 


Figure  3-3:  A  comparison  between  full  {N  =  4005)  and  reduced  outputs  (u  =  212) 
for  sensors  Si  and  S2  using  initial  condition  depicted  in  Figure  3-2  with  k  =  0.01, 
At  =  0.005  and  hundred  time  steps.  The  error  is  £  =  1.0232e“'^. 


minant  concentration  at  sensor  location  j  through  time.  The  true  initial  condition 
Uq  that  was  chosen  to  present  the  results  is  shown  in  Figure  3-5.  The  observations 
dobs  _  (^([obs^  _  _  _ ,  dg®)  were  obtained  by  propagating  uq  through  the  forward  model 
and  adding  5%  noise  to  each  component.  The  forward  model  is  given  by 

dj  =  GjUo  =  {CjA^^F)uo  for  1  <  j  <  Q,  (3.6) 

where  A  and  F  are  dehned  as  in  (2.17)  and  Cj  is  the  matrix  corresponding  to  sensor 
location  j.  As  data  uncertainties  are  Gaussian  we  can  write  the  probability  density 
distributions  of  the  values  of  contaminant  concentration  as 

Po{di,  ...,dQ)~  exp(-—  -  dff(dj  -  df)).  (3.7) 

^  i=i 
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Size  of  reduced-order  model 

Figure  3-4:  The  upper  plot  shows  the  increase  in  computational  time  in  seconds  of  the 
forward  problem  as  the  size  of  the  reduced-order  model  grows  over  the  rectangular 
domain.  The  lower  plot  depicts  the  decrease  of  the  error  dehned  in  (3.5)  as  the  size 
of  the  reduced-order  model  increases. 

From  this  we  obtain  the  probabilistic  solution  to  the  inverse  problem,  which  is  the 
posterior  probability  density  of  the  model  parameters  given  by 

(Tm{uo)  ~  exp{-^  ^{GjUo  -  df‘'Y{GjUo  -  df"")).  (3.8) 

i=i 

This  posterior  probability  density  provides  us  with  information  on  how  likely  each 
initial  scenario  uo  is.  In  the  deterministic  approach,  the  solution  of  the  inverse  problem 
is  just  a  single  estimate  of  uq,  while  the  solution  in  the  Bayesian  approach  provides 
the  probability  distribution  function  of  Uq. 

Here  we  have  two  sensors  Si  =  (0.15,0.2)  and  S2  =  (0.5,  0.2),  hence  Q  =  2.  Let 
us  analyze  this  solution  and  the  quality  of  reduced-order  models  via  Figures  3-5,  3-6, 
3-8  and  Table  3.2.  The  mean  of  the  initial  condition  held  Uq  obtained  by  the  full 
model  is  referred  to  as  and  depicted  in  the  upper  plot  of  Figure  3-6,  whereas  the 
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Figure  3-5:  This  plot  depicts  the  true  initial  contaminant  concentration  uq  used  in 
the  inverse  problem. 


mean  of  the  initial  condition  field  computed  by  the  reduced  model  of  size  n  =  212, 
referred  to  as  is  shown  in  the  lower  plot  of  Figure  3-6.  It  takes  37.6  seconds  to 
compute  using  the  full  model  and  4.2  seconds  to  compute  which  is  nine 
times  faster.  When  using  a  reduced-order  model  of  size  n  =  38  we  can  compute 
about  123  times  faster.  In  Table  3.2  and  Figure  3-8  we  present  some  reduced- 
order  model  results.  The  time  to  compute  is  referred  to  as  c  —  time,  and  the 
error  between  the  estimate  of  the  initial  condition  field  using  the  full  and  the  reduced 
model  is  denoted  by  Suo  and  defined  by 


£ 


UQ 


u, 


full 


u 


red 


(3.9) 


The  full  system  has  spatial  size  N  =  4005  and  hence  Uq  =  (ug, . . .  ,Uq,  . . .  ,Uq),  where 
Mg  refers  to  the  initial  contaminant  concentration  at  grid  point  i.  For  the  remainder 
of  this  section  i  =  2506  which  corresponds  to  the  location  (0.1477,0.2105)  in  hi.  In 
Figure  3-7  we  depict  the  marginal  probability  density  of  Ug,  namely  cr^(Mg)  which  is 
a  Gaussian  with  standard  deviation  a  and  mean  /i,  where  /x  =  =  0.2105  with  i 
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corresponding  to  grid  point  (0.1477,0.2105). 


Full  model  initial  condition  mean 
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Figure  3-6:  The  upper  plot  shows  the  full  model  estimate  of  the  initial  condition 
held.  The  lower  plot  shows  the  reduced  model  estimate  Uq^'^  of  the  initial  condition 
held  with  n  =  212  which  yields  =  0.0062. 


3.4  Reduced- Order  Model  Performance  in  Opti¬ 

mization  Problem 

3.4.1  Unconstrained  Optimization 

Let  us  consider  the  unconstrained  optimization  problem  (2.47)  introduced  in  Section 
2.4.  Due  to  the  lack  of  sensor  constraints  the  sensors  are  free  to  move  anywhere  in 
the  domain  D  dehned  in  (3.1).  The  output  region  of  interest  Qj  is  as  stated  in  (3.3) 
and  we  expect  the  optimal  sensor  locations  to  be  within  Qj.  For  the  following  results 
we  assume  the  number  of  sensors  Q  =  2.  Solving  the  optimization  problem  using 
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Values  of 


Figure  3-7:  Marginal  probability  density  with  mean  fi  =  0.2105  and  with 

varying  standard  deviation  a.  Note  that  i  corresponds  to  grid  point  (0.1477,0.2105) 
in  fl. 


the  full  model  of  size  N  =  4005  yields  optimal  locations  of  SI  =  (0.6783,  0.2034)  and 
S2  =  (0.7751,0.1999)  after  about  31.4  hours  of  computation  time.  Let  us  the  dehne 
the  average  optimization  error  Sopt  as  follows 


^  opt 


distj 


Q 


distj,  where 


+  (S'* 


(3.10) 

(3.11) 


where  optimal  sensor  locations  computed  by  the  full  model  are  S*  = 
and  optimal  sensor  locations  computed  by  the  reduced-order  model  are  given  by 
Sjj.  =  (S'*f,S'*f).  In  Table  3.3  we  present  the  optimization  results  obtained  using 
reduced-order  models  of  size  n.  The  average  optimization  error  is  dehned  in  (3.10) 
and  the  time  in  seconds  to  compute  an  optimal  solution  is  denoted  by  c  —  time. 

From  the  obtained  results  we  can  see  that  the  average  optimization  error  decreases 
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5 


Figure  3-8:  The  upper  plot  depicts  the  increase  in  computational  time  as  the  size 
of  the  reduced-order  model  increases.  The  lower  plot  depicts  the  error  according  to 
reduced-order  model  sizes  in  the  inverse  problem. 

with  increasing  n  but  only  by  very  little.  Hence  we  suggest  to  run  the  optimization 
problem  with  a  reasonably  small  reduced-order  model  as  the  error  doesn’t  change  that 
drastically  but  the  computing  time  increases  quite  fast.  Note  that  in  the  Dynamic 
Sensor  Steering  Algorithm  the  optimization  problem  takes  by  far  the  most  time  com¬ 
pared  to  the  forward  and  inverse  problem.  A  reduced-order  model  of  size  n  =  212 
computes  the  optimal  sensor  locations  about  159  times  faster  than  the  full  model  and 
with  n  =  38  we  can  compute  up  to  300  times  faster.  We  depict  the  obtained  optimal 
solutions  in  Figure  3-9. 

3.4.2  Constrained  Optimization 

In  this  section  we  consider  the  constrained  optimization  problem  (2.49)  introduced 
in  Section  2.4  with  a  radius  R  of  0.2.  The  output  region  of  interest  D/  is  as  stated  in 
(3.3)  and  we  expect  the  optimal  sensor  locations  to  be  either  within  D/,  if  the  sensor 
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case 

c  —  time 

n 

A 

1 

0.2887 

0.305 

38 

0.5 

o 

1 

2 

0.2613 

0.381 

53 

0.5 

10-® 

3 

0.0848 

0.624 

80 

0.1 

10“^ 

4 

0.0386 

1.276 

118 

0.1 

10-^^ 

5 

0.0286 

1.224 

121 

0.01 

o 

1 

6 

0.0062 

4.218 

212 

0.001 

10-^^ 

Table  3.2:  Properties  of  various  reduced-order  models  of  a  full-scale  system  with  size 
N  =  4005  and  two  output  sensors  in  the  inverse  problem. 


case 

^opt 

c  —  time 

n 

1 

0.1613 

373.921 

38 

2 

0.0338 

438.442 

53 

3 

0.0267 

651.141 

80 

4 

0.0263 

657.382 

118 

5 

0.0119 

711.401 

212 

Table  3.3:  Results  of  various  reduced-order  models  in  the  unconstrained  optimization 
problem. 


constraints  allow  it,  or  to  move  in  the  direction  of  hi/  until  the  sensor  constraint 
becomes  active.  For  the  following  results  we  assume  the  number  of  sensors  Q  =  2  with 
initial  locations  =  (0.9,  0.3)  and  S'®  =  (0.4,  0.1).  Solving  the  optimization  problem 
using  the  full  model  of  size  N  =  4005  yields  optimal  locations  of  S'l  =  (0.7549,  0.2033) 
and  5*2  =  (0.5815,  0.1840)  after  about  17  hours  of  computation  time.  In  Table  3.4  we 
present  the  optimization  results  obtained  using  reduced-order  models  of  size  n.  The 
average  optimization  error  Sopt  is  dehned  in  (3.10)  and  the  time  to  compnte  an  optimal 
solntion  is  denoted  by  c  —  time  in  seconds.  As  in  the  unconstrained  case,  the  results 
show  that  increasing  the  size  of  the  rednced-order  model  yields  only  small  reductions 
in  error,  but  has  a  large  impact  on  computing  time.  In  this  case,  a  reduced-order 
model  of  size  n  =  212  computes  the  optimal  sensor  locations  abont  27  times  faster 
than  the  full  model,  while  n  =  38  leads  to  a  speed-up  factor  of  about  94.  We  depict 
the  obtained  optimal  solutions  in  Figure  3-10. 
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Figure  3-9:  This  figure  contains  the  new  sensor  locations  computed  by  the  full  model 
and  several  reduced-order  models  in  the  unconstrained  case. 

3.4.3  Optimality  Conditions 

In  the  unconstrained  case  the  sensor  can  move  anywhere  within  the  bounds  and 
hence  we  expect  the  optimal  sensor  locations  to  be  within  Qj.  After  the  optimization 
algorithm  terminates  we  verify  that  the  gradient  at  the  new  sensor  locations  is  zero. 
Hence  we  can  state  that  we  have  found  at  least  a  local  minimum.  How  confidently  can 
we  say  that  this  local  minimum  is  in  fact  a  global  minimum.  When  we  fix  Q  — 1  sensors 
and  compute  the  objective  function,  i.e.  the  maximum  eigenvalue  of  the  covariance 
matrix,  for  the  remaining  sensor  over  the  entire  domain  hi  we  can  observe  that  the 
local  minimum  is  a  global  one  for  the  remaining  sensor.  We  can  repeat  this  procedure 
Q  times  always  letting  one  sensor  run  freely  and  we  get  similar  results.  Moreover  we 
solved  the  two  sensor  optimization  problem  for  hundred  randomly  chosen  initial  sensor 
locations  and  in  96%  of  the  cases  the  SQP  algorithm  converged  to  the  same  solution. 
In  the  remaining  4%  the  algorithm  got  stuck  near  to  the  initialized  points  because 
they  were  chosen  too  close  to  the  bonndaries  of  the  domain. 

In  the  constrained  case  we  obtain  a  solution  that  either  has  zero  gradient,  if  we 
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Figure  3-10:  This  figure  contains  the  new  sensor  locations  computed  by  the  full  model 
and  several  reduced-order  models  in  the  constrained  case. 


case 

^opt 

c  —  time 

n 

1 

0.0361 

655.644 

38 

2 

0.0302 

660.103 

53 

3 

0.0122 

688.332 

80 

4 

0.0119 

1957.057 

118 

5 

0.0086 

2304.569 

212 

Table  3.4:  Results  of  various  reduced-order  models  in  the  constrained  optimization 
problem. 


can  reach  fl/,  or  nonzero  gradients,  if  the  sensor  constraints  become  active.  We 
observed  that,  when  repeating  the  constrained  optimization  problem  long  enough, 
always  initializing  with  the  previous  optimal  solution  we  end  up  at  the  same  locations 
in  Qi  like  we  would  have  in  the  unconstrained  case.  This  is  only  true  if  all  the  sensors 
can  actually  reach  Qj.  In  the  next  section,  we  demonstrate  how  the  sensors  move 
during  various  cycles  of  the  online  stage. 
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3.5  Moving  Window  Illustration 


In  this  section  we  demonstrate  the  work  flow  of  the  Dynamic  Sensor  Steering  Algo¬ 
rithm  over  several  cycles.  We  solve  the  optimization  problem  using  a  reduced-order 
model  of  size  n  =  212.  The  quality  of  the  reduced  solution  compared  to  the  full 
solution  has  already  been  discussed  in  Section  3.4.  Note  that  we  are  also  applying  a 
set  of  sensor  constraints,  which  were  discussed  in  Section  2.4,  and  the  corresponding 
optimization  problem  is  stated  in  (2.49).  The  radius  i?  is  0.1  and  VLj  is  as  defined 
in  (3.3).  Each  figure  shows  the  current  contaminant  concentration,  the  current  sen¬ 
sor  locations,  the  circle  within  which  each  sensor  can  steer,  and  the  optimal  sensor 
locations  computed  by  solving  the  optimization  problem. 

In  Figures  3-11-3-14  we  observe  four  subsequent  cycles  of  the  online  stage  of 
the  Dynamic  Sensor  Steering  Algorithm.  Figure  3-11  corresponds  to  Cycle  1  of  the 
online  stage  and  shows  the  current  contaminant  concentration  at  tf  =  0.2,  the  do¬ 
main  of  interest  D/,  the  current  sensor  locations  =  (0.3, 0.2),  =  (0.4, 0.3), 

the  sensor  constraints  and  the  optimal  sensor  locations  S'J  =  (0.3962,0.1727)  and 
5*2  =  (0.4889,  0.2541).  The  lower  plot  in  Figure  3-11  shows  a  low-rank  approximation 
of  the  variance  held.  In  Cycle  1  the  variance  in  Qj  is  still  quite  large.  Thus  in  the 
following  cycles  we  move  the  sensors  even  closer  to  that  region.  Note  that  the  variance 
is  very  small  around  the  sensor  locations. 

Figure  3-12  corresponds  to  Cycle  2  of  the  online  stage  and  shows  the  contaminant 
concentration  at  tf  =  0.4,  the  current  sensor  locations  Aj,  computed  in  Cycle 

1,  the  sensor  constraints  and  the  optimal  sensor  locations  Sf  =  (0.4895,  0.2087)  and 
S'!  =  (0.5779,0.2093).  Due  to  the  newly  obtained  sensor  locations,  the  variance  in 
D/  has  decreased. 

Figure  3-13  corresponds  to  Cycle  3  of  the  online  stage  and  shows  the  contaminant 
concentration  at  tf  =  0.6,  the  current  sensor  locations  Sf,  Af  computed  in  Cycle 

2,  the  sensor  constraints  and  the  optimal  sensor  locations  Sf  =  (0.5895,  0.2074)  and 
S'!  =  (0.6775,0.2008).  The  variance  over  D/  decreases  further. 

Figure  3-14  corresponds  to  Cycle  4  of  the  online  stage  and  shows  the  contaminant 
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concentration  at  tj  =  0.8,  the  current  sensor  locations  Sf,  S'!  computed  in  Cycle 
3,  the  sensor  constraints  and  the  optimal  sensor  locations  Sf  =  (0.6492,0.1980)  and 
S'!  =  (0.7500,0.2002).  The  variance  held  over  Qj  approaches  zero. 

In  the  hrst  three  hgures  we  can  observe  that  the  optimal  sensor  locations  he 
on  the  boundary  imposed  by  the  sensor  constraints  and  thus  do  not  lead  to  a  zero 
gradient.  In  Figure  3-14  the  optimal  sensor  locations  both  lie  within  Qj  as  expected 
and  furthermore  none  of  the  sensor  constraints  are  active.  Moreover  with  sensors  at 
Sf  and  S'!  the  gradient  is  zero  guaranteeing  that  we  are  at  a  local  minimum. 
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Figure  3-11:  The  upper  hgure  depicts  the  contaminant  concentration  at  tj  =  0.2, 
current  and  optimal  sensor  locations  of  cycle  one  of  the  online  stage  with  constant 
velocity  held  vi.  The  lower  plot  shows  the  corresponding  variance  held  based  on  the 
current  optimal  sensor  locations  over 
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Figure  3-12:  The  upper  figure  depicts  contaminant  concentration  at  =  0.4,  current 
and  optimal  sensor  locations  of  cycle  two  of  the  online  stage  with  constant  velocity 
field  Vi.  The  lower  plot  shows  the  corresponding  variance  field  based  on  the  current 
optimal  sensor  locations  over  hi. 
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Figure  3-13:  The  upper  figure  depicts  the  contaminant  concentration  at  t/  =  0.6, 
current  and  optimal  sensor  locations  of  cycle  three  of  the  online  stage  with  constant 
velocity  field  Vi.  The  lower  plot  shows  the  corresponding  variance  field  based  on  the 
current  optimal  sensor  locations  over  fl. 
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Figure  3-14:  The  upper  figure  depicts  the  contaminant  concentration  at  tf  =  0.8, 
current  and  optimal  sensor  locations  of  cycle  four  of  the  online  stage  with  constant 
velocity  field  Vi.  The  lower  plot  shows  the  corresponding  variance  field  based  on  the 
current  optimal  sensor  locations  over  fl. 
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Chapter  4 


Model  Order  Reduction  of 
Convection-Diffusion  Equation 
with  Parameterized  Velocity  Field 


In  Chapter  4  we  consider  the  convection-diffusion  equation  and  introduce  a  parame¬ 
terized  velocity  field  in  the  Dynamic  Sensor  Steering  Algorithm  setting,  which  enables 
a  more  realistic  simulation  of  physical  processes.  The  posed  challenge  is  to  maintain 
the  algorithm’s  computational  real-time  property,  by  building  reduced-order  models 
that  are  not  only  dependent  on  the  initial  condition  but  also  on  the  velocity  held.  In 
Section  4.1  we  will  provide  a  detailed  problem  statement,  followed  by  a  description 
of  the  extended  Dynamic  Sensor  Steering  Algorithm  in  Section  4.2  and  by  a  discus¬ 
sion  of  the  applied  model  order  reduction  methodology  for  linear  dependence  of  the 
velocity  held  on  the  parameter  vector  in  Section  4.3. 

4.1  Problem  Description 

We  consider  the  the  convection-dihusion  equation  described  in  Section  2.1.2  but  now 
the  held  variable,  denoting  the  contaminant  concentration,  and  the  velocity  held,  and 
thus  the  state  solution,  depend  on  an  input  parameter  vector  /j,.  That  is,  we  now 
have  v^n)  and  ^  G  D  C  The  dihusivity  constant  is  k  and  the  observed  time 
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horizon  is  given  by  tf.  We  imposed  a  homogeneous  Dirichlet  boundary  condition  on 
Td  X  (0,t/)  and  homogeneous  Neumann  boundary  conditions  on  Fat  x  (0,t/). 

It  is  assumed  that  there  exists  a  way  of  computing  a  velocity  field  v{iJ,)  over  the 
entire  domain  hi  given  an  input  parameter  vector  fi.  Hence  there  are  two  separate 
cases,  namely  the  linear  case  if  the  dependence  of  v{iJ,)  on  fj,  is  linear  and  the  non¬ 
linear  case  when  this  dependence  is  non-linear.  In  this  chapter  we  will  focus  on  the 
linear  case. 

For  the  linear  case  we  assume  that  P  different  convective  velocity  helds  are  given 
Vi, ...  ,vp  and  that  v{iJ,)  is  an  linear  combination  of  velocity  helds  that  we  specihed 
a  priori.  Therefore  the  current  velocity  held  in  our  domain  H  can  be  computed 
by: 

p 

j=i 

with  input  parameter  /j,  =  (/ii, . . .  ,/ip)  C  M.^.  The  P  pre-computed  velocity  helds 
Vi, ...  ,Vp  could  come  from  e.g.  a  Navier-Stokes  how  solver,  or  a  reduced-order  model 
of  a  how  solver,  etc. 


4.1.1  Weak  Form  and  Finite  Element  Method 

The  weak  formulation  of  the  problem  stated  above  is  as  follows:  V/i  G  P,  hud  u(/i)  G 
X,  where  X  =  {u  G  H^{Q)  \  u  |r£,=  0},  such  that  Vtc  G  X, 

+  a{u{n),w;  fi)  =  l{w),  (4.2) 

where  a(-,  •)  is  a  bilinear  functional,  /(■)  is  a  linear  functional  and  ,w)  is 

a  bilinear  form  that  results  from  the  inner  product  of  and  the  test  functions 
w.  The  time  discretization  is  then  performed  using  Backward  Euler  or  any  other 
appropriate  scheme. 
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For  w,w  E  X  the  bilinear  and  linear  functionals  stated  in  (4.2)  are  given  by 


a(tc,  w;  ^) 

=  ai{w,w;  +  a2{w,w), 

(4.3) 

ai(tc,  ti);  fi) 

=  w  ■  Vw  dQ, 

(4.4) 

Jq 

a2{w,w) 

=  K  Xw  ■  Xw  dVt, 

(4.5) 

Jn 

l{w) 

=  [  wf  dT, 

Jq 

(4.6) 

where  k  corresponds  to  the  diffusivity  and  /  corresponds  to  the  source  term,  which 
is  zero  for  our  problem. 


Time  Discretization  and  Matrix  Eqnations 

In  Section  2.2.2  we  have  introduced  the  general  linear  discrete-time  system  (2.12)- 
(2.14),  which  can  be  expressed  in  the  following  continuous  form 


A4  u{k)  =  A  u{k),  (4.7) 

d{k)  =  C  u{k)  (4.8) 

where  it  denotes  the  vector  of  state  derivatives  with  respect  to  time  and  A4  G 
and  A  G  correspond  to  the  mass  matrix  and  the  stiffness  matrix  respectively 

obtained  from  the  hnite-element  discretization.  For  the  discretization  in  time  we  will 
use  a  Backward  Euler  method  with  a  time  step  At.  In  order  to  write  (4.7)-(4.8)  in 
matrix  form  let  us  dehne  the  following  matrices  according  to  the  chosen  Backward 
Euler  method. 


Di  - 

(4.9) 

(4.10) 

Inserting  Di  from  (4.10),  D2  from  (4.10)  and  C  from  (4.8)  into  (2.17)  we  can  derive 
the  same  matrix  form  as  stated  in  (2.15)-(2.16). 
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Streamline-Upwind  /  Petrov-Galerkin  Stabilization 


For  small  values  of  the  diffusivity  k  a  solution  to  the  discretization  of  (4.2)  can  exhibit 
non-physical  oscillations  due  to  the  negative  numerical  diffusion  of  the  Galerkin  hnite 
element  method,  which  is  described  in  detail  in  the  literature.  In  [15]  this  instability 
occurs  when  the  Peclet  number  is  greater  than  one.  In  our  problem  however  it  can  not 
be  guaranteed  to  have  the  same  Peclet  number  for  each  element  in  the  triangulation 
because  the  current  velocity  field  v  may  not  be  constant.  Thus  we  dehne  an  elemental 
Peclet  number  Pck  such  that 

Pek  =  for  k  =  1,  ,  Nelem,  (4.11) 

where  Neiem  corresponds  to  the  number  of  elements  in  our  mesh,  h  is  the  size  of  the 
element  and  Vk  denotes  the  velocity  at  the  centroid  of  element  k.  Therefore  when 
maxi<fc<jVg(g„^  Pcfc  >  1  we  will  have  add  stabilization  to  (4.2).  This  can  be  done  by 
deploying  the  Streamline-Upwind  Petrov-Galerkin  method  by  Brooks  and  Hughes  [8] , 
which  leads  to  the  weak  form 

em  p 

,w) +  a{u{p),w,p)  +  Iv{w;p)t'JZ{u{p);p)=1{w),  (4.12) 

k=i 

where  the  Streamline-Upwind  Petrov-Galerkin  stabilization  terms  are  given  by 

V{w]p)  =  v{p)-Vw  (4.13) 

^(u(Ai);M)  =  +  t^(M)  ■  Vu(Ai)  -  kVVm)-  (4.14) 

and  the  stabilization  parameter  r  can  be  computed  by 


r  = 

(4.15) 

K  = 

i  II  V  11^2  h 

2 

(4.16) 

e  = 

coth(Pe)  —  1/Pe. 

(4.17) 
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4.2  Dynamic  Sensor  Steering  Algorithm  with  Pa¬ 


rameterized  Velocity  Field 

All  introductory  definitions  from  Section  2.1.1  remain  unmodified.  It  should  be  em¬ 
phasized  though  that  contrary  to  Chapter  3,  where  m  =  (uq),  the  input  vector  is 
now  m  =  {uqjI-i),  where  uq  denotes  the  initial  condition  and  n  denotes  the  parame¬ 
ter  vector.  Moreover,  the  overall  purpose  of  the  Dynamic  Sensor  Steering  Algorithm 
doesn’t  change  but  additionally  to  the  methodology  described  in  Chapter  2  we  incor¬ 
porate  a  parameterized  velocity  held  into  the  problem.  Thus  we  take  into  account 
that  the  velocity  held  in  our  domain  D  is  subject  to  frequent  changes,  which  was  not 
captured  by  the  initial  design  of  the  algorithm  stated  in  Chapter  2.  Therefore  there 
are  signihcant  changes  in  the  algorithm  which  will  be  described  in  this  section. 

4.2.1  Offline  Stage  of  the  Dynamic  Sensor  Steering  Algo¬ 
rithm  with  Parameterized  Velocity  Field 

Figure  4-1  shows  the  offline  stage  of  the  algorithm.  The  details  of  each  step  are 
discussed  in  the  following: 

Step  Oa.  Build  reduced-order  model  gr{uo,  n)  of  physical  system  g{uo,  n)  and 
continue  working  with  reduced  system  throughout  the  Dynamic  Sensor  Steering  Al¬ 
gorithm.  Note  that  the  process  of  building  a  reduced-order  model,  that  is  now  not 
only  depending  on  an  initial  condition  uq  but  also  on  a  parameter  vector  g  &T)  G 
becomes  a  rather  complex,  challenging  task.  The  modus  operandi  differs  depending 
on  the  linearity  or  non-linearity  of  the  relationship  between  g  and  v{g).  The  linear 
approach  is  presented  in  Section  4.3. 

Step  Ob.  Place  the  Q  mobile  sensors  at  locations  S^, . . . ,  Sq  within  D  and  activate 
them.  Note  that  in  Step  1  of  the  Online  Stage  we  have  to  measure  two  values,  the 
current  contaminant  concentration  denoted  by  d*  =  {d\, . . . ,  dq),  where  i  corresponds 
to  the  cycle  we  are  in,  and  the  current  parameter  vector  g^,  which  provides  the  current 
velocity  held.  Later  on  this  issue  shall  be  addressed  and  some  solution  approaches 
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shall  be  proposed. 


Offline  Stage 


step  Oa.  Build  reduced-order  model  gr(uo,p)  of  physical  system  g(uo,p). 


Step  Ob.  Initialize  sensor  locations  Si°,...,Sq°. 


Figure  4-1:  Offline  stage  of  the  Dynamic  Sensor  Steering  Algorithm  with  parameter¬ 
ized  velocity  held. 


4.2.2  Online  Stage  of  the  Dynamic  Sensor  Steering  Algo¬ 
rithm  with  Parameterized  Velocity  Field 

Figure  4-2  shows  the  online  stage  of  the  algorithm.  The  details  of  each  step  are 
discussed  in  the  following: 

Step  1.  The  Q  sensors  create  measurements  of  the  environment  at  their  current 
locations  . . . ,  in  cycle  i.  Let  d*  =  {d\, . . . ,  (Pq)  be  the  obtained  data  on  the 
contaminant  concentration  at  the  Q  different  sensor  locations  and  f-P  be  the  obtained 
data  on  the  parameter  vector. 

Step  2.  Compute  the  current  velocity  held  u®  using  fp . 

Step  3.  Based  on  d*  =  . . . ,  dg),  pP,  u®  and  the  known  reduced  forward  model 

gr{uo,  n)  (or  Gr{uo,  g)  in  the  linear  case)  we  solve  the  inverse  problem  to  obtain  a 
prediction  of  the  initial  value  uq.  By  applying  model  order  reduction  to  g{uo,  g)  in 
the  offline  stage  we  made  sure  to  maintain  real-time  computations  throughout  the 
online  stage. 

Step  4-  In  many  applications  some  regions  Qj  within  the  observed  domain  D  need 
to  be  watched  more  thoroughly  than  others  at  any  time,  meaning  that  the  physical 
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process  that  we  want  to  predict  interests  us  especially  in  those  regions,  hi/  can  be 
hxed,  e.g.  representing  an  urban  region,  or  could  change  from  cycle  to  cycle,  e.g.  if 
we  wish  to  track  regions  of  high  concentration. 

Step  5.  We  optimize  the  sensor  locations  within  hi  such  that  the  uncertainty  in 
the  prediction  for  the  output  regions  of  interest  hi/  is  minimized.  Here,  the  full  and 
reduced  forward  models  depend  on  both  uq  and  fi. 

Step  6.  The  last  step  in  the  algorithm  steers  the  mobile  sensors  from  their  current 
locations  . . . ,  S^^^  to  the  optimized  locations  SI, Sq  within  Tgteer  and  return 
to  Step  1. 

4.3  Model  Order  Reduction  Methodology  in  the 
Linear  Case 

In  this  section  we  assume  that  the  dependence  of  the  convective  velocity  held  v{pL)  on 
the  parameter  vector  ^  is  linear.  Now  let  us  observe  once  more  the  weak  formulation 
derived  in  Section  4.1.1  and  pay  particular  attention  to  the  functional  n)  which 
depends  on  the  parameter  vector  /j,.  In  equation  (4.2)  we  divided  functional  a(-,  ■;  n) 
into  two  independent  functionals  ai(-,-;/i)  and  a2(-,-)  where  the  latter  corresponds 
to  the  diffusive  part  of  the  problem  and  ai(-,  ■;  ^)  shows  the  dependence  on  fi  that 
is  inherent  to  the  convective  part  of  the  problem.  Due  to  their  independence  of  fi, 
both  /(■)  and  a2(-,  •)  can  be  computed  offline  once,  leading  to  the  mass  matrix  A4  and 
the  right-hand  side  vector  JF  that  will  not  change  throughout  the  problem  (assuming 
constant  diffusivity  k).  Note  that  T  is  zero  in  our  problem.  Conversely, 

ai{w,w]  i-i)  =  /  w  v{pl)  ■  Vw  dVL  \/  w,w  E  X  (4-18) 

Jn 

shows  that  we  cannot  compute  ai(-,  ■;  /x)  once,  because  /j,  and  thus  the  stiffness  matrix 
A  are  subject  to  frequent  changes.  Potentially,  this  could  mean  that  we  can  no  longer 
compute  a  reduced-order  model  in  the  offline  Stage  of  the  Dynamic  Sensor  Steering 
Algorithm.  The  consequence  is  a  loss  of  executing  the  algorithm  in  real-time.  However 
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let  us  now  incorporate  (4.1)  into  (4.18),  which  yields 


w  v{fi)  ■  Ww  dVt  = 


/  ^  d-j  ^j)  ■  Vh;  did 

,=i 

^  r 

/ 

,=i 


w  Vj  ■  Vtc  dil. 


(4.19) 

(4.20) 


From  this  we  can  see  that  the  integral  over  hi  in  (4.20)  is  no  longer  dependent  on 
and  can  therefore  be  computed  offline  for  each  Vj,  j  =  1, . . . ,  P.  This  means  that  we 
now  have  to  pre-compute  P  different  stiffness  matrices  . . .  ,Ap  corresponding  to 
the  P  different  velocity  helds. 


Thus,  the  continuous  representation  of  the  discrete-time  system  stated  in  (4.7)- 
(4.8)  becomes 

M  iiik)  =  A{ii)  u(k),  (4.21) 

d{k)  =  C  u{k),  (4.22) 

p 

with  ^(^)  =  Aj,  (4.23) 

i=i 

and  from  the  Backward  Euler  method  for  time  discretization  we  obtain 

1  ^ 

Pf  =  (4.24) 

D;  =  (4.25) 

If  we  dehne  the  matrices 


I 

0 

0 

c 

0 

...  0 

-D2 

0 

0 

c 

0 

A^^  = 

0 

-P2 

pf 

,  C  = 

0 

c 

0 

0 

0 

0 

-P2  Pf 

0 

0  c 

68 


(4.26) 


and 


I  r  1  r 

m(0)  d{0) 

m(1)  d{l) 

0  ,  u=  ,  d= 

u(T)  d(T) 

oj 

the  new  matrix  form  of  the  system  becomes: 

u  =  F  Mo,  (4.28) 

d  =  C  u.  (4.29) 

Since  we  have  decoupled  /j,  from  the  model  matrices,  we  can  compute  all  necessary 
reduced-order  model  matrices  offline.  The  reduced  matrices  are  dehned  by 


Mr  =  V^MV 

(4.30) 

Aj^r  =  V'^AjV  j  =  1, . . . ,  P 

(4.31) 

Cr  =  CV, 

(4.32) 

where  computation  of  the  basis  V  was  discussed  in  Section  2.2.2. 

model,  (4.24)-(4.25)  becomes 

In  the  reduced 

^  p 

i=i 

(4.33) 

(4.34) 

which  ultimately  leads  to  the  reduced  system  in  matrix  form 

Af:*  Ur  =  Fr  Mo, 

(4.35) 

dr  r  ^iJjr , 

(4.36) 

(4.27) 
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Online  Stage 
Cycle  i 


Figure  4-2:  One  cycle  of  the  online  stage  of  the  Dynamic  Sensor  Steering  Algorithm 
with  parameterized  velocity  held. 
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Chapter  5 


Application  to 

Convection-Diffusion  Systems  with 
Parameterized  Velocity  Fields 


In  Chapter  5  we  present  the  results  obtained  by  the  Dynamic  Sensor  Steering  Algo¬ 
rithm  with  parameterized  velocity  held  in  the  linear  case  described  in  Section  4.3. 
First,  let  us  compare  reduced-order  model  and  full  model  performance  over  a  rec¬ 
tangular  domain  in  Section  5.1.  Then  in  Section  5.2,  we  apply  the  algorithm  to  a 
backwards-facing  step  domain  with  velocity  helds  computed  by  a  Navier-Stokes  how 
solver. 


5.1  Dynamic  Sensor  Steering  over  Rectangular  Do¬ 
main 

5.1.1  Description  of  Rectangular  Domain 

The  rectangular  domain  used  throughout  this  section  was  described  in  Section  3.1. 
Note  that  boundary  conditions  as  well  as  the  output  region  of  interest  remain  un¬ 
changed  and  that  all  following  results  are  obtained  assuming  two  sensors.  As  previ¬ 
ously  mentioned  in  Section  4.1.1  we  need  to  stabilize  our  system  in  case  the  maximum 
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elemental  Peclet  number  Pe^,  defined  in  (4.11),  exceeds  one.  Hence,  let  us  choose  a 
diffusivity  k  ensuring 

max  Pcfc  <  1  (5.1) 

by  finding  an  upper  bound  for  (4.11).  We  consider  only  low  Peclet  number  flows  in 
order  to  avoid  stabilization  because  otherwise  we  would  have  to  pre-compute  stabiliza¬ 
tion  matrices  which  are  described  in  Section  4.1.1  and  introduce  additional  notation 
here. 


5.1.2  A  Priori  Velocity  Fields  over  Rectangular  Domain 


The  Dynamic  Sensor  Steering  Algorithm  in  the  affine  case  requires  P  pre-computed 
velocity  fields  and  an  input  parameter  ^  G  In  this  section  the  a  priori  velocity 
fields  shall  be  represented  by  polynomial  functions  in  x  and  y.  In  particular  we  want 
to  be  able  to  represent  velocity  fields  modeled  by  all  polynomial  functions  fi, . . . ,  fu, 
where  U  is  the  number  of  velocity  fields,  up  to  a  degree  of  y  which  causes 

degifu)  <X  for  M  =  1, . . . ,  P.  (5.2) 

If  there  is  no  redundancy  within  fi, . . . ,  fu  then  U  =  P.  From  this  we  immediately 
get  that  P  >  X  +  1  and  the  polynomials  of  degree  smaller  or  equal  to  y  have  a 
resulting  vector  space  of  dimension  x  -|-  1.  In  order  to  span  this  space  let  us  assume 
a  monomial  basis 


fu{x) 

=  ^  for  M  =  1, . 

•  • )  X  +  1) 

(5.3) 

fu{y) 

=  y^~^  for  M  =  1, . 

•  •,x  +  i- 

(5.4) 

74 

Then  P  has  to  be  2  *  (x  +  1)  and  let  us  define  the  corresponding  velocity  fields  as 

fkix) 

Vk  =  for  A;  =  1, . . .  ,x  +  1  (5.5) 

0 

0  1 

Vk+x+i  =  for  /c  =  1, . . . ,  X  +  1  (5.6) 

-fk{y) 

which  makes  sure  that  we  can  adequately  model  polynomial  flow  up  to  an  order  of  x 
in  both  the  x  and  y  direction. 

For  the  computational  results  over  the  rectangular  domain  we  choose  x  =  3,  which 
leads  to  P  =  8  a  priori  velocity  fields 


and  fj,  will  be  specified  at  the  beginning  of  each  new  cycle  of  the  Dynamic  Sensor 
Steering  Algorithm  and  shall  be  stated  accordingly. 

5.1.3  Moving  Window  Illustration 

As  we  have  already  discussed  the  performance  of  reduced-order  models  in  the  forward 
problem,  inverse  problem  and  optimization  problem  over  the  rectangular  domain  in 
Chapter  3  the  reader  is  referred  to  Sections  3.2  -  3.4.  Let  us  now  depict  the  work 
flow  of  the  algorithm  as  the  velocity  field  changes.  Each  figure  below  corresponds  to 
a  different  cycle  of  the  online  stage  of  the  Dynamic  Sensor  Steering  Algorithm.  The 
current  velocity  field  that  is  represented  by  the  vector  changes  after  each  cycle. 
Again  we  apply  sensor  constraints  stated  in  (2.49).  The  radius  R  is  0.1.  We  assumed 
K  =  0.01,  At  =  0.005,  unchanged  Qj  and  the  initial  condition  displayed  in  3-2. 

Figure  5-1  corresponds  to  Cycle  1  of  the  online  stage  and  shows  the  current  con- 
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taminant  concentration  dX  tj  =  0.2,  the  domain  of  interest  the  current  sensor 
locations  5'°  =  (0.7,0.35),  =  (0.4,  0.2),  the  sensor  constraints  and  the  optimal 

sensor  locations  Si  =  (0.7189,0.2518)  and  =  (0.4992,0.2123). 
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^  Current  sensor  location 
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Figure  5-1:  Contaminant  concentration  at  =  0.2,  current  and  optimal  sen¬ 
sor  locations  of  cycle  one  of  the  online  stage  with  parameterized  velocity  held  by 
H  =  {0.75,0.25,0,0,0,0,0,0}. 


Figure  5-2  corresponds  to  Cycle  2  of  the  online  stage  and  shows  the  current  con¬ 
taminant  concentration  at  =  0.4,  the  domain  of  interest  Qj,  the  current  sensor 
locations  5}  and  5}  as  computed  in  Cycle  1,  the  sensor  constraints  and  the  optimal 
sensor  locations  Sf  =  (0.7392,0.2029)  and  =  (0.5989,0.2203). 

Figure  5-3  corresponds  to  Cycle  3  of  the  online  stage  and  shows  the  current  con¬ 
taminant  concentration  at  =  0.6,  the  domain  of  interest  Qi,  the  current  sensor 
locations  Sf  and  S'!  as  computed  in  Cycle  2,  the  sensor  constraints  and  the  optimal 
sensor  locations  Sf  =  (0.7610,0.1839)  and  S'!  =  (0.6684,0.2011)  which  is  a  local 
minimum. 
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Figure  5-2:  Contaminant  concentration  ai  tf  =  0.4,  current  and  optimal  sen¬ 
sor  locations  of  cycle  two  of  the  online  stage  with  parameterized  velocity  held  by 
^  =  {0.25,  0.10, 0.25,  0,  0.15,  0,  0,  0.25}. 

5.2  Dynamic  Sensor  Steering  over  Backward  Fac¬ 
ing  Step 

5.2.1  Description  of  Backward  Facing  Step  Domain 

The  dimension  of  the  backward  facing  step  domain  is  best  described  by  (5.9)  and 
the  corresponding  computational  domain  has  spatial  mesh  size  of  iV  =  2417.  Figure 
5-4  depicts  the  subdivision  of  in  triangles.  We  have  k  =  0.01,  At  =  0.0025  and 
the  output  region  of  interest  is  Qj.  The  following  results  were  computed  assuming 
two  sensors.  At  the  inhow  boundary  x  =  0  and  0  <  y  <  0.4  we  apply  homogeneous 
Dirichlet  conditions  and  everywhere  else  homogeneous  Neumann  boundary  conditions 
are  applied. 


0.2  <y<0A  if  0<  a;  <0.2 
0<y  <0.4  if  0.2  <  a;  <  1 


(5.9) 
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Figure  5-3:  Contaminant  concentration  ai  tf  =  0.6,  current  and  optimal  sensor 
locations  of  cycle  three  of  the  online  stage  with  parameterized  velocity  held  by 
^  =  {0,0,0.25,0,0.75,0,0,0}. 


5.2.2  A  Priori  Velocity  Fields  over  Backward  Facing  Step 
Domain 


More  realistic  velocity  helds  for  the  backward  facing  step  domain  were  computed  using 
a  Navier-Stokes  how  solver,  i.e.  the  Incompressible  Flow  Iterative  Solution  Software 
(IFISS)  [16,  17],  The  Navier-Stokes  equations  can  be  written  as 

-u'V'^v  +  v-Vv  +  Vp  =  /,  (5.10) 

V-u  =  0,  (5.11) 

where  u  corresponds  to  the  kinematic  viscosity  constant,  f  corresponds  to  a  source 
term,  ,  v  is  the  velocity  held  and  p  represents  pressure.  Boundary  conditions  are 
given  by 


V  =  m  on  F^) , 


dv 

u— - np  = 

on 


N, 
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0  on  F 


(5.12) 

(5.13) 


x-axis 


Figure  5-4:  Triangular  mesh  of  the  backward  facing  step  domain.  Note  that  the  mesh 
used  for  computations  is  much  hner. 


where  n  is  the  normal  vector,  F £>  corresponds  to  the  boundary  where  we  apply  ho¬ 
mogeneous  Dirichlet  conditions  and  F at  corresponds  to  the  boundary  where  we  apply 
homogeneous  Neumann  conditions.  In  particular,  a  Poiseuille  flow  is  imposed  on  the 
inflow  boundary,  a  no-flow  (zero  velocity)  condition  is  imposed  on  the  top  and  bottom 
boundaries  and  a  Neumann  condition  is  applied  at  the  outflow  boundary  automati¬ 
cally  setting  the  mean  outflow  pressure  to  zero.  To  solve  the  Navier-Stokes  equations 
stated  in  (5.12)-(5.13)  a  hnite  element  method  using  a  quadrilateral  element  mesh  is 
applied.  The  resulting  non-linear  algebraic  system  is  solved  using  iterative  methods, 
i.e.  a  hybrid  method  performing  a  small  number  of  Picard  iterations  first  for  getting 
a  good  starting  point  for  Newton’s  method.  For  further  details  on  Krylov  subspace 
methods,  stabilization  and  pre-conditioning  please  refer  to  [16]. 

To  demonstrate  the  computational  results  over  the  backward  facing  step  domain 
we  choose  P  =  4  obtaining  velocity  helds  Ui, . . .  ,up  corresponding  to  four  different 
choices  of  the  kinematic  viscosities  specihed  in  (5.14).  Note  that  the  boundary  condi¬ 
tions  for  the  Navier-Stokes  equations  remain  unchanged  when  computing  the  solution 
corresponding  to  different  viscosities. 


z/i  =  1/50,  1/2  =  1/75,  1/3  =  1/25,  z/4  =  1/10  (5.14) 
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As  the  resulting  velocity  fields  can  not  be  described  analytically  as  in  Section  5.1.2  let 
us  present  the  solution  of  the  Navier-Stokes  equations  graphically  for  vi.  Figure  5-5 
contains  the  solution  of  the  velocity  field  by  plotting  the  streamlines  and  the  solution 
of  the  pressure  field.  In  Figure  5-6  we  show  the  velocity  at  certain  nodes  in  the  step 
domain.  Note  that  will  be  specified  at  the  beginning  of  each  new  cycle  of  the 
Dynamic  Sensor  Steering  Algorithm  but  due  to  the  fact  that  we  can  not  interpolate 
between  velocity  fields  with  different  viscosities  only  one  fij  will  be  set  to  one  whereas 
all  the  other  ones  must  remain  zero. 


Pressure  field 


Figure  5-5:  The  upper  plot  depicts  the  streamlines  over  the  backward  facing  step 
domain  computed  using  a  viscosity  z/  =  1/50.  The  lower  plot  shows  the  pressure  field 
over  the  backward  facing  step  domain  computed  using  a  viscosity  u  =  1/50. 
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Figure  5-6:  Velocity  at  selected  grid  points  of  the  backward  facing  step  domain  com¬ 
puted  using  a  viscosity  v  =  1/50. 

5.2.3  Reduced- Order  Model  Performance  in  Forward  Prob¬ 
lem 

In  this  section  we  assumed  a  diffusivity  k  =  0.01,  At  =  0.0025  and  an  input  vector 
=  {1,0,  0,0}  specifying  the  velocity  held 

=  (5.15) 

where  corresponds  to  the  viscosity  vi  dehned  in  (5.14).  We  choose  to  demonstrate 
the  results  using  a  superposition  of  hve  Gaussian  functions  as  the  initial  contaminant 
concentration  depicted  in  Figure  5-7.  Each  Gaussian  function  is  dehned  as  in  (3.4). 
In  Figure  5-8  we  compare  reduced  model  outputs  to  full-scale  outputs  at  two  sensor 
locations  Si  =  (0.15,0.25)  and  S'2  =  (0.45,0.2).  Table  5.1  displays  various  reduced- 
order  model  properties,  where  e  denotes  the  error  compared  to  the  full  model  dehned 
in  (3.5),  c  —  time  denotes  the  computational  time  to  obtain  the  outputs  in  seconds, 
n  denotes  the  size  of  the  reduced-order  model,  and  the  signihcance  of  A  and  Ji  is 
discussed  in  Section  3.2.  Running  the  forward  problem  using  the  full  model  takes 
4.88  seconds.  From  Table  5.1  we  see  that  the  largest  reduced-order  model  of  size 
n  =  185  is  about  twelve  times  faster  than  the  full  model  and  the  smallest  reduced- 
order  model  of  size  n  =  32  computes  results  even  350  times  faster.  The  performance 
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with  respect  to  error  and  computational  time  of  each  reduced-order  model  can  be 
viewed  in  Figure  5-9. 


Contaminant  Concentration  at  t=0 
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0.3  # 
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Figure  5-7:  Sample  test  initial  condition  used  to  compare  reduced  model  to  full  model 
in  the  forward  problem  over  backward  facing  step  domain. 

5.2.4  Solution  of  the  Inverse  Problem 

As  theoretical  details  have  already  been  discussed  in  Section  3.3,  we  will  only  present 
the  numerical  results  over  the  backward  facing  step  domain  here.  The  true  initial 
condition  uq  that  was  chosen  to  obtain  the  observations  dohs  is  depicted  in  Figure  5-7. 

The  mean  of  the  initial  condition  held  Uq  obtained  by  the  full  model  is  referred  to 
as  and  depicted  in  the  upper  plot  of  Figure  5-10,  whereas  the  mean  of  the  initial 
condition  held  computed  by  the  reduced  model  of  size  n  =  185,  referred  to  as 
is  shown  in  the  lower  plot  of  Figure  5-10.  It  takes  34.62  seconds  to  compute 
using  the  full  model  and  5.55  seconds  to  compute  When  using  a  reduced-order 
model  of  size  32  we  can  compute  about  150  times  faster  than  the  full  model 
computations.  In  Table  5.2  we  present  some  reduced-order  model  results.  The  time 
to  compute  is  referred  to  as  c  —  time,  and  the  error  between  the  estimate  of 
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Figure  5-8:  A  comparison  between  full  {N  =  2417)  and  reduced  outputs  (u  =  185) 
for  sensors  Si  and  S2  using  initial  condition  depicted  in  Figure  5-7  with  k  =  0.01, 

At  =  0.0025  and  100  time  steps.  The  error  is  £  =  8.5419e“^. 

the  initial  condition  held  using  the  full  and  the  reduced  model  is  denoted  by  e^o  and 
dehned  by  (3.9). 

5.2.5  Reduced-Order  Model  Performance  in  Optimization  Prob¬ 
lem 

Again  we  consider  the  unconstrained  optimization  problem  (2.47)  introduced  in  Sec¬ 
tion  2.4  with  a  velocity  held  corresponding  to  u  =  1/50  and  At  =  0.0025.  Due  to 
the  lack  of  senor  constraints  the  sensors  are  free  to  move  anywhere  in  the  domain  D 
dehned  in  (5.9).  The  output  region  of  interest  flj  is  as  stated  in  (3.3)  and  we  expect 
the  optimal  sensor  locations  to  be  within  Qj.  For  the  following  presented  results  we 
assume  the  number  of  sensors  Q  =  2.  Solving  the  optimization  problem  using  the 
full  model  of  size  N  =  2417  yields  optimal  locations  of  SI  =  (0.6463,0.1901)  and 
S2  =  (0.7100,0.1998)  after  about  19  hours.  The  SQP  algorithm  terminated  because 
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Figure  5-9:  The  upper  plot  shows  the  increase  in  computational  time  in  seconds  of 
the  forward  problem  as  the  size  of  the  reduced-order  model  grows  over  the  backward 
facing  step  domain.  The  lower  plot  depicts  the  decrease  of  the  error  defined  in  (3.5) 
as  the  size  of  the  reduced-order  model  increases. 

the  number  of  maximum  function  evaluations  was  exceeded.  The  gradient  at  the 
optimal  solution  was  reasonably  small  to  account  for  at  least  a  local  minimum.  We 
have  observed  that  when  velocity  helds  are  computed  by  a  Navier-Stokes  solver,  the 
objective  function  is  less  well  behaved.  Compared  to  results  in  Chapter  3  only  73% 
of  random  initial  locations  converge  to  the  above  optimum.  The  remaining  27%  ei¬ 
ther  get  stuck  in  local  minima  that  are  not  within  Qj  (and  yield  a  bigger  objective 
value  than  the  optimal  point)  or  when  initialized  too  close  to  a  boundary  do  not 
move  at  all  (we  assume  that  this  is  due  to  flatness  in  the  objective  function).  In 
Table  5.3  we  present  the  optimization  results  obtained  using  reduced-order  models 
of  size  n.  The  average  optimization  error  Sopt  is  defined  in  (3.10)  and  the  time  to 
compute  an  optimal  solution  is  denoted  by  c  —  time.  The  following  conclusions  can 
be  drawn  from  Table  5.3:  the  smallest  reduced-order  model  does  not  perform  the 
fastest.  Even  though  a  single  function  evaluation  using  the  reduced-order  model  of 
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case 

e 

c  —  time 

n 

A 

/i 

1 

0.0858 

0.0140 

32 

0.5 

o 

1 

2 

0.0599 

0.0210 

45 

0.5 

10-6 

3 

0.0055 

0.0570 

71 

0.1 

o 

1 

4 

0.0026 

0.1060 

105 

0.1 

10-6 

5 

0.0017 

0.0880 

96 

0.01 

10-^ 

6 

0.0016 

0.2680 

146 

0.01 

10-6 

7 

8.5419e-^ 

0.4030 

185 

0.001 

10-6 

Table  5.1:  Properties  of  various  reduced-order  models  of  a  full-scale  system  with 
size  N  =  2417  and  two  output  sensors  over  the  backward  facing  step  domain  in  the 
forward  problem.  The  error  e  is  dehned  in  (3.5)  and  the  initial  condition  used  is 
depicted  in  Fignre  5-7. 


case 

£uo 

c  —  time 

n 

A 

/i 

1 

0.7492 

0.224 

32 

0.5 

10-^ 

2 

0.7125 

0.320 

45 

0.5 

10-6 

3 

0.6070 

0.664 

71 

0.1 

10-^ 

4 

0.5264 

1.533 

105 

0.1 

10-6 

5 

0.4256 

3.000 

146 

0.01 

10-6 

5 

0.4025 

5.551 

185 

0.001 

10-6 

Table  5.2:  Properties  of  various  reduced-order  models  of  a  full-scale  system  with  size 
N  =  2417  and  two  output  sensors  in  the  inverse  problem. 


size  n  =  45  is  indeed  faster  than  using  a  rednced-order  model  of  size  n  =  105,  the 
number  of  iterations  needed  in  the  SQP  algorithm  is  not  dependent  on  the  size  of 
the  reduced-order  model.  Hence  the  fastest  result  is  not  necessarily  obtained  by  the 
smallest  reduced-order  model.  However,  in  general  (when  comparing  average  run 
times)  smaller  rednced-order  models  yield  a  smaller  compntational  time.  Moreover 
we  can  observe  that  increasing  the  size  of  the  reduced-order  model  by  four  leads  to 
an  error  that  is  rednced  by  a  factor  of  three.  Hence  we  snggest  that  using  a  smaller 
reduced-order  model  in  the  optimization  problem  yields  a  small  computational  time 
with  a  reasonably  small  error.  We  omit  details  about  the  constrained  optimization 
problem  as  we  observe  similar  behavior  as  in  the  unconstrained  case. 
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Figure  5-10:  The  upper  plot  shows  the  full  model  estimate  of  the  initial  condition 
held  over  the  backward  facing  step  domain.  The  lower  plot  shows  the  reduced  model 
estimate  of  the  initial  condition  held  with  n  =  185  which  yields  =  0.4025. 


5.2.6  Moving  Window  Illustration 

In  this  section  we  demonstrate  the  work  how  over  the  backward  facing  step  domain. 
We  solve  the  optimization  problem  using  a  reduced-order  model  of  size  n  =  185  and 
apply  a  set  of  sensor  constraints  stated  in  (2.49).  The  radius  R  that  was  used  here  is 
0.15. 

Figure  5-11  corresponds  to  Cycle  1  of  the  online  stage  and  shows  the  current 
contaminant  concentration  at  tj  =  0.2,  the  domain  of  interest  Qj,  the  current  sensor 
locations  =  (0.3, 0.2),  =  (0.9, 0.1),  the  sensor  constraints  and  the  optimal 

sensor  locations  S'J  =  (0.4436,0.2434)  and  S'g  =  (0.7714,0.1772).  We  also  include  a 
plot  corresponding  to  the  not  normalized  variance  held.  In  Cycle  1  the  variance  on 
the  right  boundary  of  Qj  is  very  small  because  one  sensor  is  located  there.  On  the 
left  boundary  however  the  variance  is  still  big  as  the  other  sensor  is  still  to  far  away 


case 

^opt 

c  —  time 

n 

1 

0.0295 

262.061 

45 

2 

0.0195 

242.479 

105 

3 

0.0166 

278.449 

146 

4 

0.0105 

464.317 

185 

Table  5.3:  Results  of  various  reduced-order  models  in  the  unconstrained  optimization 
problem  over  the  backward-facing  step. 

from  Qj. 

Figure  5-12  corresponds  to  Cycle  2  of  the  online  stage  and  shows  the  current 
contaminant  concentration  at  tf  =  0.4,  the  domain  of  interest  Qj,  the  current  sensor 
locations  S'!  and  as  computed  in  Cycle  1,  the  sensor  constraints  and  the  optimal 
sensor  locations  Sf  =  (0.5887,0.2814)  and  S'!  =  (0.6992,0.1987).  The  variance  in  Qj 
is  close  to  zero  now  that  one  sensor  is  in  the  output  region  of  interest  and  the  other 
one  is  very  close  to  it. 

Figure  5-13  corresponds  to  Cycle  3  of  the  online  stage  and  shows  the  current 
contaminant  concentration  at  tf  =  0.6,  the  domain  of  interest  Qj,  the  current  sensor 
locations  Sf  and  S'!  as  computed  in  Cycle  2,  the  sensor  constraints  and  the  optimal 
sensor  locations  Sf  =  (0.7113,0.1950)  and  Rf  =  (0.6511,0.1740).  As  we  decrease  the 
variance  in  Qj  by  moving  both  sensors  into  this  region,  the  variance  at  regions  that 
are  further  away  from  the  sensor  locations  increases.  Once  again  let  us  mention  that 
we  only  want  to  minimize  uncertainty  in  predictions  over  Qj  and  hence  neglect  the 
rest  of  the  domain. 
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Figure  5-11:  The  upper  figure  depicts  the  contaminant  concentration  at  tf  =  0.2,  cur¬ 
rent  and  optimal  sensor  locations  of  cycle  one  of  the  online  stage  with  parameterized 
velocity  field  hy  =  {1,0,  0,0}.  The  lower  figure  shows  the  corresponding  variance 
field  based  on  the  current  optimal  sensor  locations  over  fl. 
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Figure  5-12:  The  upper  figure  depicts  the  contaminant  concentration  at  t/  =  0.4,  cur¬ 
rent  and  optimal  sensor  locations  of  cycle  two  of  the  online  stage  with  parameterized 
velocity  field  by  /^  =  {0,  0, 1,  0}.  The  lower  figure  shows  the  corresponding  variance 
field  based  on  the  current  optimal  sensor  locations  over  fl. 
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Figure  5-13:  The  upper  figure  depicts  the  contaminant  concentration  at  tf  =  0.6, 
current  and  optimal  sensor  locations  of  cycle  three  of  the  online  stage  with  parame¬ 
terized  velocity  field  by  ^  =  {0,0,0,!}.  The  lower  figure  shows  the  corresponding 
variance  field  based  on  the  current  optimal  sensor  locations  over  fl. 
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Chapter  6 


Conclusions  and  Future  Work 

6.1  Summary  and  Conclusions 

In  this  thesis  we  develop  the  methodology  for  a  real-time  Dynamic  Sensor  Steering 
Algorithm  that  computes  the  best-possible  prediction  of  an  ongoing  physical  process 
by  combining  model  order  reduction  with  a  Bayesian  approach  to  inverse  problems 
and  with  optimization.  We  apply  the  algorithm  to  a  contaminant  transport  problem, 
paying  particular  attention  to  full  model  versus  reduced-order  model  performance. 

The  algorithm  is  divided  into  two  separate  stages:  the  offline  stage,  which  is  exe¬ 
cuted  only  once,  and  the  online  stage,  which  consists  of  a  measure-predict-optimize- 
steer  cycle  and  is  executed  repeatedly.  In  the  offline  stage  we  place  mobile  sensors 
in  the  domain  that  we  wish  to  observe  and  we  build  a  reduced-order  model  of  the 
physical  system,  using  the  reduced  basis  technique,  which  projects  the  large-scale  sys¬ 
tem  onto  a  basis  in  a  space  of  reduced  dimensions.  This  basis  is  obtained  by  proper 
orthogonal  decomposition  and  Hessian-based  model  reduction. 

In  the  online  stage  we  first  obtain  observations  of  the  contaminant  concentration 
at  current  sensor  locations  by  measurements  and  then  we  solve  an  inverse  problem 
yielding  a  prediction  of  the  initial  contamination  using  a  Bayesian  approach.  The 
solution  to  the  inverse  problem  is  a  probability  density  function  of  the  initial  contam¬ 
inant  concentration.  This  probability  density  function  provides  us  with  information 
on  how  likely  each  possible  contamination  scenario  is  based  on  the  obtained  mea- 
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surements.  The  knowledge  of  the  contamination  distribution  provides  us  with  the 
advantage  of  being  able  to  instruct  hre  hghters  and  police  more  effectively  and  if 
necessary  even  evacuate  people  who  are  in  danger.  The  physical  process,  i.e.  the 
contamination,  changes  with  time  and  hence  we  move  the  sensors  to  new  locations 
that  minimize  the  uncertainty  in  the  prediction.  We  obtain  those  new  locations  by 
solving  an  optimization  problem,  with  non-linear  constraints  if  sensor  constraints  are 
imposed,  using  the  sequential  quadratic  programming  algorithm.  After  the  sensors 
are  steered  to  their  new  locations,  the  online  stage  repeats  itself. 


In  reality  the  velocity  held,  e.g.  the  wind  velocity,  is  also  subject  to  changes,  which 
alters  the  way  the  contamination  evolves.  Hence  we  extend  the  Dynamic  Sensor 
Steering  Algorithm  to  handle  a  parameterized  velocity  held,  i.e.  we  pre-compute 
several  diherent  velocity  helds,  using  a  Navier-Stokes  how  solver  and  obtain  diherent 
physical  systems  according  to  diherent  velocity  helds.  In  the  offline  stage  we  build 
reduced-order  models  for  each  of  the  physical  systems.  In  the  online  stage  we  are  now 
capable  of  approximating  the  full  model  by  using  the  reduced-model,  with  a  velocity 
held  that  is  any  linear  combination  of  the  previously  computed  velocity  helds,  enabling 
to  model  a  much  broader  range  of  scenarios  accurately.  The  performance  in  terms 
of  computational  time  and  accuracy  of  reduced-order  models  versus  full  models  in 
the  forward  problem,  the  inverse  problem  and  the  optimization  problem  is  discussed 
individually.  Due  to  reduced-order  modeling  we  are  able  to  execute  one  cycle  of 
the  online  stage  about  260  times  faster  than  with  the  full  model  with  average  relative 
errors  of  magnitude  0(10“^).  The  methodology  applied  to  the  contaminant  transport 
problem  is  widely  applicable  to  all  sorts  of  problems  where  we  wish  to  observe  certain 
phenomena  whose  location  or  features  are  not  known  before  hand.  The  following 
section  includes  some  suggestions  for  improvements  and  future  work. 
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6.2  Future  Work 


6.2.1  Recommendations 

In  this  thesis  the  general  methodology  was  applied  to  a  linear  two-dimensional  con¬ 
taminant  transport  problem.  It  is  a  challenging  and  interesting  task  to  extend  the 
methodology  to  three-dimensional  problems  and  non-linear  problems. 

In  case  the  dependence  of  the  velocity  held  on  the  input  parameter  /j,  in  the 
Dynamic  Sensor  Steering  Algorithm  with  parameterized  velocity  held  is  nonlinear,  the 
velocity  held  could  be  anything  but  in  order  to  reduce  the  level  of  abstraction  we 
can  think  of  this  model  as  a  Navier-Stokes  how  solver  that  computes  the  velocity  held 
for  every  point  in  the  domain  D  given  fi.  Then  ^  could  contain  a  specihc  Reynolds 
number,  boundary  conditions,  initial  conditions,  etc.  We  can  extend  the  algorithm  to 
additionally  sample  for  obtaining  current  velocity  helds  in  each  cycle  instead  of  pre- 
computing  velocity  helds.  An  approach  for  this  extension  is  proposed  in  the  following 
section. 

In  references  [34,  35]  a  technique  based  on  proper  orthogonal  decomposition  is 
proposed  that  reconstructs  a  close  approximation  to  the  velocity  held.  In  this  thesis 
we  have  focused  on  computing  reduced-order  models  to  accurately  approximate  the 
full  model  of  a  contaminant  transport  problem  while  assuming  certain  velocity  helds. 
Combining  the  work  that  was  presented  in  [34,  35]  and  this  thesis  can  yield  even  more 
efficient  and  more  realistic  results. 

6.2.2  Model  Order  Reduction  Methodology  in  the  Nonlinear 
Case 

We  consider  the  more  general  case  where  the  dependence  of  parameter  vector  on  the 
convective  velocity  held  v{^)  is  non-linear.  In  order  to  avoid  online  re-computation 
of  the  reduced-order  model  we  need  to  obtain  an  affine  decomposition  of  the  bilin¬ 
ear  functional  ai{w,w;  n)  in  terms  of  products  of  parameter-dependent  coefficients, 
for  1  <  m  <  n  (computed  online),  and  parameter-independent  bilinear  forms 
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a'^{w,w)  for  1  <  m  <  n  (computed  offline).  The  resulting  bilinear  functional  can  be 
written  as 

n 

=  (6.1) 

m=l 

In  order  to  achieve  this  we  introduce  a  parameter  sample  that  adequately  covers  the 
parameter  space  V, 

Sk  =  . . . ,  hk  e  T>}.  (6.2) 

Then  we  compute  solution  snapshots 

for  1  <  A;  <  K}.  (6.3) 

Then  we  apply  POD  to  the  snapshots  and  hence  obtain  the  basis  functions 
for  our  reduced-order  model  and  dehne  the  approximation  space  V  = 

span{(i,. . .  ,Cn}-  For  the  following  discussion  let  us  assume  the  basis  is  orthonor- 
malized.  Now  we  also  generate  snapshots  of  the  velocity  held, 

‘^K  =  {n(A^fc)  for  1  <  /c  <  K}.  (6.4) 

Then  we  use  the  ’’Best  Points”  Interpolation  Method  (BPIM)  introduced  in  [40]  to 
construct  a  set  of  interpolation  points  the  basis  functions  such 

that  we  can  dehne  a  coefflcient-function  approximation,  for  the  convective 

velocity  held  v{x;  /j,)  as  follows, 

n 

Vn{x-  =  v(z^;  ii)  (6.5) 

m=l 

where  x  =  {x,y)  G  D.  Then  with  ipjizi)  =  Sij,  where  6ij  is  the  Kronecker 

symbol,  can  thus  be  obtained  by 

n 

Cii^)  =  l<i<n.  (6.6) 

m=l 


94 


Then  the  discrete  version  of  the  bilinear  functional  ai{w,w,  n)  is  given  by 


dn+^V2{zl;fi)  f  dVl 

m=i  dn  Jn  y 

(67) 

n  n 

=  <x(0,Ci)  +  a^(Ci,Ci)  i<bi<«,  (6.8) 

m=l  m=l 

where  v{z'^]  /j,)  =  (ni(z^;  ^),  m))-  Now  we  can  solve  the  weak  form  (4.2)  in  the 

reduced  basis  space  V.  The  same  affine  decomposition  described  above  for  ai(tc,  w;  /j,) 
should  also  be  used  for  all  the  SUPG  stabilization  terms  in  (4.2). 

When  actually  implementing  this  method  we  need  to  address  the  following  issue. 
Throughout  the  algorithm  the  sensors  that  measure  the  ongoing  contamination  steer 
to  locations  which  minimize  uncertainty  in  the  prediction,  however  the  functionality 
of  the  best  interpolation  points  that  we  use  to  approximate  ai(-,  ■;  ^)  is  completely 
unrelated  to  the  sensor  locations.  Thus  we  can  either  introduce  a  second  set  of  sensors 
that  merely  depends  on  the  best  interpolation  points  or  we  can  only  work  with  one  set 
of  sensors  and  try  to  hnd  a  reasonable  balance  between  the  locations  that  minimize 
uncertainty  and  locations  of  best  interpolation  points.  Although  it  is  more  costly  to 
deploy  another  set  of  sensors  the  results  will  to  be  more  accurate.  However  it  should 
be  considered  that  when  working  with  one  set  of  sensors  due  to  limited  resources 
hnding  the  optimal  trade-off  between  minimizing  uncertainty  in  the  prediction  and 
approximating  the  functional  ai(-,-;/i)  is  a  very  interesting  problem  by  itself  and 
should  be  examined  in  future  work. 
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